25class PseudoDifferentialOperator:
  26    """
  27    Pseudo-differential operator with dynamic symbol evaluation on spatial grids.
  28    Supports both 1D and 2D operators, and can be defined explicitly (symbol mode)
  29    or extracted automatically from symbolic equations (auto mode).
  30
  31    Parameters
  32    ----------
  33    expr : sympy expression
  34        Symbolic expression representing the pseudo-differential symbol.
  35    vars_x : list of sympy symbols
  36        Spatial variables (e.g., [x] for 1D, [x, y] for 2D).
  37    var_u : sympy function, optional
  38        Function u(x, t) used in auto mode to extract the operator symbol.
  39    mode : str, {'symbol', 'auto'}
  40        - 'symbol': directly uses expr as the operator symbol.
  41        - 'auto': computes the symbol automatically by applying expr to exp(i x ξ).
  42
  43    Attributes
  44    ----------
  45    dim : int
  46        Spatial dimension (1 or 2).
  47    fft, ifft : callable
  48        Fast Fourier transform and inverse (scipy.fft or scipy.fft2).
  49    p_func : callable
  50        Evaluated symbol function ready for numerical use.
  51
  52    Notes
  53    -----
  54    - In 'symbol' mode, `expr` should be expressed in terms of spatial variables and frequency variables (ξ, η).
  55    - In 'auto' mode, the symbol is derived by applying the differential expression to a complex exponential.
  56    - Frequency variables are internally named 'xi' and 'eta' for consistency.
  57    - Uses numpy for numerical evaluation and scipy.fft for FFT operations.
  58
  59    Examples
  60    --------
  61    >>> # Example 1: 1D Laplacian operator (symbol mode)
  62    >>> from sympy import symbols
  63    >>> x, xi = symbols('x xi', real=True)
  64    >>> op = PseudoDifferentialOperator(expr=xi**2, vars_x=[x], mode='symbol')
  65
  66    >>> # Example 2: 1D transport operator (auto mode)
  67    >>> from sympy import Function
  68    >>> u = Function('u')
  69    >>> expr = u(x).diff(x)
  70    >>> op = PseudoDifferentialOperator(expr=expr, vars_x=[x], var_u=u(x), mode='auto')
  71    """
  72
  73    def __init__(self, expr, vars_x, var_u=None, mode='symbol'):
  74        self.dim = len(vars_x)
  75        self.mode = mode
  76        self.symbol_cached = None
  77        self.expr = expr
  78        self.vars_x = vars_x
  79
  80        if self.dim == 1:
  81            x, = vars_x
  82            xi_internal = symbols('xi', real=True)
  83            expr = expr.subs(symbols('xi', real=True), xi_internal)
  84            self.fft = partial(fft, workers=FFT_WORKERS)
  85            self.ifft = partial(ifft, workers=FFT_WORKERS)
  86
  87            if mode == 'symbol':
  88                self.p_func = lambdify((x, xi_internal), expr, 'numpy')
  89                self.symbol = expr
  90            elif mode == 'auto':
  91                if var_u is None:
  92                    raise ValueError("var_u must be provided in mode='auto'")
  93                exp_i = exp(I * x * xi_internal)
  94                P_ei = expr.subs(var_u, exp_i)
  95                symbol = simplify(P_ei / exp_i)
  96                symbol = expand(symbol)
  97                self.symbol = symbol
  98                self.p_func = lambdify((x, xi_internal), symbol, 'numpy')
  99            else:
 100                raise ValueError("mode must be 'auto' or 'symbol'")
 101
 102        elif self.dim == 2:
 103            x, y = vars_x
 104            xi_internal, eta_internal = symbols('xi eta', real=True)
 105            expr = expr.subs(symbols('xi', real=True), xi_internal)
 106            expr = expr.subs(symbols('eta', real=True), eta_internal)
 107            self.fft = partial(fft2, workers=FFT_WORKERS)
 108            self.ifft = partial(ifft2, workers=FFT_WORKERS)
 109
 110            if mode == 'symbol':
 111                self.symbol = expr
 112                self.p_func = lambdify((x, y, xi_internal, eta_internal), expr, 'numpy')
 113            elif mode == 'auto':
 114                if var_u is None:
 115                    raise ValueError("var_u must be provided in mode='auto'")
 116                exp_i = exp(I * (x * xi_internal + y * eta_internal))
 117                P_ei = expr.subs(var_u, exp_i)
 118                symbol = simplify(P_ei / exp_i)
 119                symbol = expand(symbol)
 120                self.symbol = symbol
 121                self.p_func = lambdify((x, y, xi_internal, eta_internal), symbol, 'numpy')
 122            else:
 123                raise ValueError("mode must be 'auto' or 'symbol'")
 124
 125        else:
 126            raise NotImplementedError("Only 1D and 2D supported")
 127
 128        if mode == 'auto':
 129            print("\nsymbol = ")
 130            pprint(self.symbol, num_columns=NUM_COLS)
 131        
 132    def evaluate(self, X, Y, KX, KY, cache=True):
 133        """
 134        Evaluate the pseudo-differential operator's symbol on a grid of spatial and frequency coordinates.
 135
 136        The method dynamically selects between 1D and 2D evaluation based on the spatial dimension.
 137        If caching is enabled and a cached symbol exists, it returns the cached result to avoid recomputation.
 138
 139        Parameters
 140        ----------
 141        X, Y : ndarray
 142            Spatial grid coordinates. In 1D, Y is ignored.
 143        KX, KY : ndarray
 144            Frequency grid coordinates. In 1D, KY is ignored.
 145        cache : bool, default=True
 146            If True, stores the computed symbol for reuse in subsequent calls to avoid redundant computation.
 147
 148        Returns
 149        -------
 150        ndarray
 151            Evaluated symbol values over the input grid. Shape matches the input spatial/frequency grids.
 152
 153        Raises
 154        ------
 155        NotImplementedError
 156            If the spatial dimension is not 1D or 2D.
 157        """
 158        if cache and self.symbol_cached is not None:
 159            return self.symbol_cached
 160
 161        if self.dim == 1:
 162            symbol = self.p_func(X, KX)
 163        elif self.dim == 2:
 164            symbol = self.p_func(X, Y, KX, KY)
 165
 166        if cache:
 167            self.symbol_cached = symbol
 168
 169        return symbol
 170
 171    def clear_cache(self):
 172        """
 173        Clear cached symbol evaluations.
 174        """        
 175        self.symbol_cached = None
 176
 177    def apply(self, u, x_grid, kx, boundary_condition='periodic', 
 178              y_grid=None, ky=None, dealiasing_mask=None,
 179              freq_window='gaussian', clamp=1e6, space_window=False):
 180        """
 181        Apply the pseudo-differential operator to the input field u.
 182    
 183        This method dispatches the application of the pseudo-differential operator based on:
 184        
 185        - Whether the symbol is spatially dependent (x/y)
 186        - The boundary condition in use (periodic or dirichlet)
 187    
 188        Supported operations:
 189        
 190        - Constant-coefficient symbols: applied via Fourier multiplication.
 191        - Spatially varying symbols: applied via Kohn–Nirenberg quantization.
 192        - Dirichlet boundary conditions: handled with non-periodic convolution-like quantization.
 193    
 194        Dispatch Logic:\n
 195        if not self.is_spatial: u ↦ Op(p)(D) ⋅ u = 𝓕⁻¹[ p(ξ) ⋅ 𝓕(u) ]\n
 196        elif periodic: u ↦ Op(p)(x,D) ⋅ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ based of FFT (quicker)\n
 197        elif dirichlet: u ↦ Op(p)(x,D) ⋅ u ≈ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ (slower)\n
 198        
 199        Parameters
 200        ----------
 201        u : ndarray
 202            Function to which the operator is applied
 203        x_grid : ndarray
 204            Spatial grid in x direction
 205        kx : ndarray
 206            Frequency grid in x direction
 207        boundary_condition : str
 208            'periodic' or 'dirichlet'
 209        y_grid : ndarray, optional
 210            Spatial grid in y direction (for 2D)
 211        ky : ndarray, optional
 212            Frequency grid in y direction (for 2D)
 213        dealiasing_mask : ndarray, optional
 214            Dealiasing mask
 215        freq_window : str
 216            Frequency windowing ('gaussian' or 'hann')
 217        clamp : float
 218            Clamp symbol values to [-clamp, clamp]
 219        space_window : bool
 220            Apply spatial windowing
 221            
 222        Returns
 223        -------
 224        ndarray
 225            Result of applying the operator
 226        """
 227        # Check if symbol depends on spatial variables
 228        is_spatial = self._is_spatial_dependent()
 229        
 230        # Case 1: Constant symbol with periodic BC (fast path)
 231        if not is_spatial and boundary_condition == 'periodic':
 232            return self._apply_constant_fft(u, x_grid, kx, y_grid, ky, dealiasing_mask)
 233        
 234        # Case 2: Spatial symbol with periodic BC
 235        elif boundary_condition == 'periodic':
 236            symbol_func = self._get_symbol_func()
 237            return kohn_nirenberg_fft(
 238                u_vals=u,
 239                symbol_func=symbol_func,
 240                x_grid=x_grid,
 241                kx=kx,
 242                fft_func=self.fft,
 243                ifft_func=self.ifft,
 244                dim=self.dim,
 245                y_grid=y_grid,
 246                ky=ky,
 247                freq_window=freq_window,
 248                clamp=clamp,
 249                space_window=space_window
 250            )
 251        
 252        # Case 3: Dirichlet BC (non-periodic)
 253        elif boundary_condition == 'dirichlet':
 254            symbol_func = self._get_symbol_func()
 255            
 256            if self.dim == 1:
 257                return kohn_nirenberg_nonperiodic(
 258                    u_vals=u,
 259                    x_grid=x_grid,
 260                    xi_grid=kx,
 261                    symbol_func=symbol_func,
 262                    freq_window=freq_window,
 263                    clamp=clamp,
 264                    space_window=space_window
 265                )
 266            elif self.dim == 2:
 267                return kohn_nirenberg_nonperiodic(
 268                    u_vals=u,
 269                    x_grid=(x_grid, y_grid),
 270                    xi_grid=(kx, ky),
 271                    symbol_func=symbol_func,
 272                    freq_window=freq_window,
 273                    clamp=clamp,
 274                    space_window=space_window
 275                )
 276        
 277        else:
 278            raise ValueError(f"Invalid boundary condition '{boundary_condition}'")
 279    
 280    def _is_spatial_dependent(self):
 281        """
 282        Check if the symbol depends on spatial variables.
 283        
 284        Returns
 285        -------
 286        bool
 287            True if symbol depends on x (or x, y)
 288        """
 289        if self.dim == 1:
 290            return self.symbol.has(self.vars_x[0])
 291        elif self.dim == 2:
 292            x, y = self.vars_x
 293            return self.symbol.has(x) or self.symbol.has(y)
 294        else:
 295            return False
 296    
 297    def _get_symbol_func(self):
 298        """
 299        Get a lambdified version of the symbol.
 300        
 301        Returns
 302        -------
 303        callable
 304            Lambdified symbol function
 305        """
 306        if self.dim == 1:
 307            x = self.vars_x[0]
 308            xi = symbols('xi', real=True)
 309            return lambdify((x, xi), self.symbol, 'numpy')
 310        elif self.dim == 2:
 311            x, y = self.vars_x
 312            xi, eta = symbols('xi eta', real=True)
 313            return lambdify((x, y, xi, eta), self.symbol, 'numpy')
 314        else:
 315            raise NotImplementedError("Only 1D and 2D supported")
 316    
 317    def _apply_constant_fft(self, u, x_grid, kx, y_grid, ky, dealiasing_mask):
 318        """
 319        Apply a constant-coefficient pseudo-differential operator in Fourier space.
 320
 321        This method assumes the symbol is diagonal in the Fourier basis and acts as a 
 322        multiplication operator. It performs the operation:
 323        
 324            (ψu)(x) = 𝓕⁻¹[ -σ(k) · 𝓕[u](k) ]
 325
 326        where:
 327        - σ(k) is the combined pseudo-differential operator symbol
 328        - 𝓕 denotes the forward Fourier transform
 329        - 𝓕⁻¹ denotes the inverse Fourier transform
 330
 331        The dealiasing mask is applied before returning to physical space.
 332        
 333        Parameters
 334        ----------
 335        u : ndarray
 336            Input function
 337        x_grid : ndarray
 338            Spatial grid (x)
 339        kx : ndarray
 340            Frequency grid (x)
 341        y_grid : ndarray, optional
 342            Spatial grid (y, for 2D)
 343        ky : ndarray, optional
 344            Frequency grid (y, for 2D)
 345        dealiasing_mask : ndarray, optional
 346            Dealiasing mask
 347            
 348        Returns
 349        -------
 350        ndarray
 351            Result
 352        """
 353        u_hat = self.fft(u)
 354        
 355        # Evaluate symbol at grid points
 356        if self.dim == 1:
 357            X_dummy = np.zeros_like(kx)
 358            symbol_vals = self.p_func(X_dummy, kx)
 359        elif self.dim == 2:
 360            KX, KY = np.meshgrid(kx, ky, indexing='ij')
 361            X_dummy = np.zeros_like(KX)
 362            Y_dummy = np.zeros_like(KY)
 363            symbol_vals = self.p_func(X_dummy, Y_dummy, KX, KY)
 364        else:
 365            raise ValueError("Only 1D and 2D supported")
 366        
 367        # Apply symbol
 368        u_hat *= symbol_vals
 369        
 370        # Apply dealiasing
 371        if dealiasing_mask is not None:
 372            u_hat *= dealiasing_mask
 373        
 374        return self.ifft(u_hat)
 375
 376    def principal_symbol(self, order=1):
 377        """
 378        Compute the leading homogeneous component of the pseudo-differential symbol.
 379
 380        This method extracts the principal part of the symbol, which is the dominant 
 381        term under high-frequency asymptotics (|ξ| → ∞). The expansion is performed 
 382        in polar coordinates for 2D symbols to maintain rotational symmetry, then 
 383        converted back to Cartesian form.
 384
 385        Parameters
 386        ----------
 387        order : int
 388            Order of the asymptotic expansion in powers of 1/ρ, where ρ = |ξ| in 1D 
 389            or ρ = sqrt(ξ² + η²) in 2D. Only the leading-order term is returned.
 390
 391        Returns
 392        -------
 393        sympy.Expr
 394            The principal symbol component, homogeneous of degree `m - order`, where 
 395            `m` is the original symbol's order.
 396
 397        Notes:
 398        - In 1D, uses direct series expansion in ξ.
 399        - In 2D, expands in radial variable ρ while preserving angular dependence.
 400        - Useful for microlocal analysis and constructing parametrices.
 401        """
 402
 403        p = self.symbol
 404        if self.dim == 1:
 405            xi = symbols('xi', real=True, positive=True)
 406            return simplify(series(p, xi, oo, n=order).removeO())
 407        elif self.dim == 2:
 408            xi, eta = symbols('xi eta', real=True, positive=True)
 409            # Homogeneous radial expansion: we set (ξ, η) = ρ (cosθ, sinθ)
 410            rho, theta = symbols('rho theta', real=True, positive=True)
 411            p_rho = p.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 412            expansion = series(p_rho, rho, oo, n=order).removeO()
 413            # Revert back to (ξ, η)
 414            expansion_cart = expansion.subs({rho: sqrt(xi**2 + eta**2),
 415                                             cos(theta): xi / sqrt(xi**2 + eta**2),
 416                                             sin(theta): eta / sqrt(xi**2 + eta**2)})
 417            return simplify(powdenest(expansion_cart, force=True))
 418                       
 419    def is_homogeneous(self, tol=1e-10):
 420        """
 421        Check whether the symbol is homogeneous in the frequency variables.
 422    
 423        Returns
 424        -------
 425        (bool, Rational or float or None)
 426            Tuple (is_homogeneous, degree) where:
 427            - is_homogeneous: True if the symbol satisfies p(λξ, λη) = λ^m * p(ξ, η)
 428            - degree: the detected degree m if homogeneous, or None
 429        """
 430        from sympy import symbols, simplify, expand, Eq
 431        from sympy.abc import l
 432    
 433        if self.dim == 1:
 434            xi = symbols('xi', real=True, positive=True)
 435            l = symbols('l', real=True, positive=True)
 436            p = self.symbol
 437            p_scaled = p.subs(xi, l * xi)
 438            ratio = simplify(p_scaled / p)
 439            if ratio.has(xi):
 440                return False, None
 441            try:
 442                deg = simplify(ratio).as_base_exp()[1]
 443                return True, deg
 444            except Exception:
 445                return False, None
 446    
 447        elif self.dim == 2:
 448            xi, eta = symbols('xi eta', real=True, positive=True)
 449            l = symbols('l', real=True, positive=True)
 450            p = self.symbol
 451            p_scaled = p.subs({xi: l * xi, eta: l * eta})
 452            ratio = simplify(p_scaled / p)
 453            # If ratio == l**m with no (xi, eta) left, it's homogeneous
 454            if ratio.has(xi, eta):
 455                return False, None
 456            try:
 457                base, exp = ratio.as_base_exp()
 458                if base == l:
 459                    return True, exp
 460            except Exception:
 461                pass
 462            return False, None
 463
 464    def symbol_order(self, max_order=10, tol=1e-3):
 465        """
 466        Estimate the homogeneity order of the pseudo-differential symbol in high-frequency asymptotics.
 467    
 468        This method attempts to determine the leading-order behavior of the symbol p(x, ξ) or p(x, y, ξ, η)
 469        as |ξ| → ∞ (in 1D) or |(ξ, η)| → ∞ (in 2D). The returned value represents the asymptotic growth or decay rate,
 470        which is essential for understanding the regularity and mapping properties of the corresponding operator.
 471    
 472        The function uses symbolic preprocessing to ensure proper factorization of frequency variables,
 473        especially in sqrt and power expressions, to avoid erroneous order detection (e.g., due to hidden scaling).
 474    
 475        Parameters
 476        ----------
 477        max_order : int, optional
 478            Maximum number of terms to consider in the series expansion. Default is 10.
 479        tol : float, optional
 480            Tolerance threshold for evaluating the coefficient magnitude. If the coefficient is too small,
 481            the detected order may be discarded. Default is 1e-3.
 482    
 483        Returns
 484        -------
 485        float or None
 486            - If the symbol is homogeneous, returns its exact homogeneity degree as a float.
 487            - Otherwise, estimates the dominant asymptotic order from leading terms in the expansion.
 488            - Returns None if no valid order could be determined.
 489    
 490        Notes
 491        -----
 492        - In 1D:
 493            Two strategies are used:
 494                1. Expand directly in xi at infinity.
 495                2. Substitute xi = 1/z and expand around z = 0.
 496    
 497        - In 2D:
 498            - Transform the symbol into polar coordinates: (xi, eta) = rho*(cos(theta), sin(theta)).
 499            - Expand in rho at infinity, then extract the leading term's power.
 500            - An alternative substitution using 1/z is also tried if the first method fails.
 501    
 502        - Preprocessing steps:
 503            - Sqrt expressions involving frequencies are rewritten to isolate the leading variable.
 504            - Power expressions are factored explicitly to ensure correct symbolic scaling.
 505    
 506        - If the symbol is not homogeneous, a warning is issued, and the result should be interpreted with care.
 507        
 508        - For non-homogeneous symbols, only the principal asymptotic term is considered.
 509    
 510        Raises
 511        ------
 512        NotImplementedError
 513            If the spatial dimension is neither 1 nor 2.
 514        """
 515        from sympy import (
 516            symbols, series, simplify, sqrt, cos, sin, oo, powdenest, radsimp,
 517            expand, expand_power_base
 518        )
 519    
 520        def preprocess_sqrt(expr, freq):
 521            return expr.replace(
 522                lambda e: e.func == sqrt and freq in e.free_symbols,
 523                lambda e: freq * sqrt(1 + (e.args[0] - freq**2) / freq**2)
 524            )
 525    
 526        def preprocess_power(expr, freq):
 527            return expr.replace(
 528                lambda e: e.is_Pow and freq in e.free_symbols,
 529                lambda e: freq**e.exp * (1 + e.base / freq**e.base.as_powers_dict().get(freq, 0))**e.exp
 530            )
 531    
 532        def validate_order(power, coeff, vars_x, tol):
 533            if power is None:
 534                return None
 535            if any(v in coeff.free_symbols for v in vars_x):
 536                print("⚠️ Coefficient depends on spatial variables; ignoring")
 537                return None
 538            try:
 539                coeff_val = abs(float(coeff.evalf()))
 540                if coeff_val < tol:
 541                    print(f"⚠️ Coefficient too small ({coeff_val:.2e} < {tol})")
 542                    return None
 543            except Exception as e:
 544                print(f"⚠️ Coefficient evaluation failed: {e}")
 545                return None
 546            return int(power) if power == int(power) else float(power)
 547    
 548        # Homogeneity check
 549        is_homog, degree = self.is_homogeneous()
 550        if is_homog:
 551            return float(degree)
 552        else:
 553            print("⚠️ The symbol is not homogeneous. The asymptotic order is not well defined.")
 554    
 555        if self.dim == 1:
 556            x = self.vars_x[0]
 557            xi = symbols('xi', real=True, positive=True)
 558    
 559            try:
 560                print("1D symbol_order - method 1")
 561                expr = preprocess_sqrt(self.symbol, xi)
 562                s = series(expr, xi, oo, n=max_order).removeO()
 563                lead = simplify(powdenest(s.as_leading_term(xi), force=True))
 564                power = lead.as_powers_dict().get(xi, None)
 565                coeff = lead / xi**power if power is not None else 0
 566                print("lead =", lead)
 567                print("power =", power)
 568                print("coeff =", coeff)
 569                order = validate_order(power, coeff, [x], tol)
 570                if order is not None:
 571                    return order
 572            except Exception:
 573                pass
 574    
 575            try:
 576                print("1D symbol_order - method 2")
 577                z = symbols('z', real=True, positive=True)
 578                expr_z = preprocess_sqrt(self.symbol.subs(xi, 1/z), 1/z)
 579                s = series(expr_z, z, 0, n=max_order).removeO()
 580                lead = simplify(powdenest(s.as_leading_term(z), force=True))
 581                power = lead.as_powers_dict().get(z, None)
 582                coeff = lead / z**power if power is not None else 0
 583                print("lead =", lead)
 584                print("power =", power)
 585                print("coeff =", coeff)
 586                order = validate_order(power, coeff, [x], tol)
 587                if order is not None:
 588                    return -order
 589            except Exception as e:
 590                print(f"⚠️ fallback z failed: {e}")
 591            return None
 592    
 593        elif self.dim == 2:
 594            x, y = self.vars_x
 595            xi, eta = symbols('xi eta', real=True, positive=True)
 596            rho, theta = symbols('rho theta', real=True, positive=True)
 597    
 598            try:
 599                print("2D symbol_order - method 1")
 600                p_rho = self.symbol.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 601                p_rho = preprocess_power(preprocess_sqrt(p_rho, rho), rho)
 602                s = series(simplify(p_rho), rho, oo, n=max_order).removeO()
 603                lead = radsimp(simplify(powdenest(s.as_leading_term(rho), force=True)))
 604                power = lead.as_powers_dict().get(rho, None)
 605                coeff = lead / rho**power if power is not None else 0
 606                print("lead =", lead)
 607                print("power =", power)
 608                print("coeff =", coeff)
 609                order = validate_order(power, coeff, [x, y], tol)
 610                if order is not None:
 611                    return order
 612            except Exception as e:
 613                print(f"⚠️ polar expansion failed: {e}")
 614    
 615            try:
 616                print("2D symbol_order - method 2")
 617                z = symbols('z', real=True, positive=True)
 618                xi_eta = {xi: (1/z) * cos(theta), eta: (1/z) * sin(theta)}
 619                p_rho = preprocess_sqrt(self.symbol.subs(xi_eta), 1/z)
 620                s = series(simplify(p_rho), z, 0, n=max_order).removeO()
 621                lead = radsimp(simplify(powdenest(s.as_leading_term(z), force=True)))
 622                power = lead.as_powers_dict().get(z, None)
 623                coeff = lead / z**power if power is not None else 0
 624                print("lead =", lead)
 625                print("power =", power)
 626                print("coeff =", coeff)
 627                order = validate_order(power, coeff, [x, y], tol)
 628                if order is not None:
 629                    return -order
 630            except Exception as e:
 631                print(f"⚠️ fallback z (2D) failed: {e}")
 632            return None
 633    
 634        else:
 635            raise NotImplementedError("Only 1D and 2D supported.")
 636
 637    
 638    def asymptotic_expansion(self, order=3):
 639        """
 640        Compute the asymptotic expansion of the symbol as |ξ| → ∞ (high-frequency regime).
 641    
 642        This method expands the pseudo-differential symbol in inverse powers of the 
 643        frequency variable(s), either in 1D or 2D. It handles both polynomial and 
 644        exponential symbols by performing a series expansion in 1/|ξ| up to the specified order.
 645    
 646        The expansion is performed directly in Cartesian coordinates for 1D symbols.
 647        For 2D symbols, the method uses polar coordinates (ρ, θ) to perform the expansion 
 648        at infinity in ρ, then converts the result back to Cartesian coordinates.
 649    
 650        Parameters
 651        ----------
 652        order : int, optional
 653            Maximum order of the asymptotic expansion. Default is 3.
 654    
 655        Returns
 656        -------
 657        sympy.Expr
 658            The asymptotic expansion of the symbol up to the given order, expressed in Cartesian coordinates.
 659            If expansion fails, returns the original unexpanded symbol.
 660    
 661        Notes:
 662        - In 1D: expansion is performed directly in terms of ξ.
 663        - In 2D: the symbol is first rewritten in polar coordinates (ρ,θ), expanded asymptotically 
 664          in ρ → ∞, then converted back to Cartesian coordinates (ξ,η).
 665        - Handles special case when the symbol is an exponential function by expanding its argument.
 666        - Symbolic normalization is applied early (via `simplify`) for 2D expressions to improve convergence.
 667        - Robust to failures: catches exceptions and issues warnings instead of raising errors.
 668        - Final expression is simplified using `powdenest` and `expand` for improved readability.
 669        """
 670        p = self.symbol
 671    
 672        if self.dim == 1:
 673            xi = symbols('xi', real=True, positive=True)
 674    
 675            try:
 676                # Case: exponential function
 677                if p.func == exp and len(p.args) == 1:
 678                    arg = p.args[0]
 679                    arg_series = series(arg, xi, oo, n=order).removeO()
 680                    expanded = series(exp(expand(arg_series)), xi, oo, n=order).removeO()
 681                    return simplify(powdenest(expanded, force=True))
 682                else:
 683                    expanded = series(p, xi, oo, n=order).removeO()
 684                    return simplify(powdenest(expanded, force=True))
 685    
 686            except Exception as e:
 687                print(f"Warning: 1D expansion failed: {e}")
 688                return p
 689    
 690        elif self.dim == 2:
 691            xi, eta = symbols('xi eta', real=True, positive=True)
 692            rho, theta = symbols('rho theta', real=True, positive=True)
 693    
 694            # Normalize before substitution
 695            p = simplify(p)
 696    
 697            # Substitute polar coordinates
 698            p_polar = p.subs({
 699                xi: rho * cos(theta),
 700                eta: rho * sin(theta)
 701            })
 702    
 703            try:
 704                # Handle exponentials
 705                if p_polar.func == exp and len(p_polar.args) == 1:
 706                    arg = p_polar.args[0]
 707                    arg_series = series(arg, rho, oo, n=order).removeO()
 708                    expanded = series(exp(expand(arg_series)), rho, oo, n=order).removeO()
 709                else:
 710                    expanded = series(p_polar, rho, oo, n=order).removeO()
 711    
 712                # Convert back to Cartesian
 713                norm = sqrt(xi**2 + eta**2)
 714                expansion_cart = expanded.subs({
 715                    rho: norm,
 716                    cos(theta): xi / norm,
 717                    sin(theta): eta / norm
 718                })
 719    
 720                # Final simplifications
 721                result = simplify(powdenest(expansion_cart, force=True))
 722                result = expand(result)
 723                return result
 724    
 725            except Exception as e:
 726                print(f"Warning: 2D expansion failed: {e}")
 727                return p  
 728            
 729    def compose_asymptotic(self, other, order=1, mode='kn', sign_convention=None):
 730        """
 731        Compose two pseudo-differential operators using an asymptotic expansion
 732        in the chosen quantization scheme (Kohn–Nirenberg or Weyl).
 733    
 734        Parameters
 735        ----------
 736        other : PseudoDifferentialOperator
 737            The operator to compose with this one.
 738        order : int, default=1
 739            Maximum order of the asymptotic expansion.
 740        mode : {'kn', 'weyl'}, default='kn'
 741            Quantization mode:
 742            - 'kn' : Kohn–Nirenberg quantization (left-quantized)
 743            - 'weyl' : Weyl symmetric quantization
 744        sign_convention : {'standard', 'inverse'}, optional
 745            Controls the phase factor convention for the KN case:
 746            - 'standard' → (i)^(-n), gives [x, ξ] = +i (physics convention)
 747            - 'inverse' → (i)^(+n), gives [x, ξ] = -i (mathematical adjoint convention)
 748            If None, defaults to 'standard'.
 749    
 750        Returns
 751        -------
 752        sympy.Expr
 753            Symbolic expression for the composed symbol up to the given order.
 754    
 755        Notes
 756        -----
 757        - In 1D (Kohn–Nirenberg):
 758            (p ∘ q)(x, ξ) ~ Σₙ (1/n!) (i sgn)^n ∂_ξⁿ p(x, ξ) ∂_xⁿ q(x, ξ)
 759        - In 1D (Weyl):
 760            (p # q)(x, ξ) = exp[(i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q)] p(x, ξ) q(x, ξ)
 761            truncated at given order.
 762    
 763        Examples
 764        --------
 765        X = a*x, Y = b*ξ
 766        X_op.compose_asymptotic(Y_op, order=3, mode='weyl')
 767        """
 768    
 769        from sympy import diff, factorial, simplify, symbols
 770    
 771        assert self.dim == other.dim, "Operator dimensions must match"
 772        p, q = self.symbol, other.symbol
 773    
 774        # Default sign convention
 775        if sign_convention is None:
 776            sign_convention = 'standard'
 777        sign = -1 if sign_convention == 'standard' else +1
 778    
 779        # --- 1D case ---
 780        if self.dim == 1:
 781            x = self.vars_x[0]
 782            xi = symbols('xi', real=True)
 783            result = 0
 784    
 785            if mode == 'kn':  # Kohn–Nirenberg
 786                for n in range(order + 1):
 787                    term = (1 / factorial(n)) * diff(p, xi, n) * diff(q, x, n) * (1j) ** (sign * n)
 788                    result += term
 789    
 790            elif mode == 'weyl':  # Weyl symmetric composition
 791                # Weyl star product: exp((i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q))
 792                result = 0
 793                for n in range(order + 1):
 794                    for k in range(n + 1):
 795                        # k derivatives acting as (∂_ξ^k p)(∂_x^(n−k) q)
 796                        coeff = (1 / (factorial(k) * factorial(n - k))) * ((1j / 2) ** n) * ((-1) ** (n - k))
 797                        term = coeff * diff(p, xi, k, x, n - k, evaluate=True) * diff(q, x, k, xi, n - k, evaluate=True)
 798                        result += term
 799    
 800            else:
 801                raise ValueError("mode must be either 'kn' or 'weyl'")
 802    
 803            return simplify(result)
 804    
 805        # --- 2D case ---
 806        elif self.dim == 2:
 807            x, y = self.vars_x
 808            xi, eta = symbols('xi eta', real=True)
 809            result = 0
 810    
 811            if mode == 'kn':
 812                for n in range(order + 1):
 813                    for i in range(n + 1):
 814                        j = n - i
 815                        term = (1 / (factorial(i) * factorial(j))) * \
 816                               diff(p, xi, i, eta, j) * diff(q, x, i, y, j) * (1j) ** (sign * n)
 817                        result += term
 818    
 819            elif mode == 'weyl':
 820                for n in range(order + 1):
 821                    for i in range(n + 1):
 822                        j = n - i
 823                        coeff = (1 / (factorial(i) * factorial(j))) * ((1j / 2) ** n) * ((-1) ** (n - i))
 824                        term = coeff * diff(p, xi, i, eta, j, x, 0, y, 0) * diff(q, x, i, y, j, xi, 0, eta, 0)
 825                        result += term
 826            else:
 827                raise ValueError("mode must be either 'kn' or 'weyl'")
 828    
 829            return simplify(result)
 830    
 831        else:
 832            raise NotImplementedError("Only 1D and 2D cases are implemented")
 833
 834    def commutator_symbolic(self, other, order=1, mode='kn', sign_convention=None):
 835        """
 836        Compute the symbolic commutator [A, B] = A∘B − B∘A of two pseudo-differential operators
 837        using formal asymptotic expansion of their composition symbols.
 838    
 839        This method computes the asymptotic expansion of the commutator's symbol up to a given 
 840        order, based on the symbolic calculus of pseudo-differential operators in the 
 841        Kohn–Nirenberg quantization. The result is a purely symbolic sympy expression that 
 842        captures the leading-order noncommutativity of the operators.
 843    
 844        Parameters
 845        ----------
 846        other : PseudoDifferentialOperator
 847            The pseudo-differential operator B to commute with this operator A.
 848        order : int, default=1
 849            Maximum order of the asymptotic expansion. 
 850            - order=1 yields the leading term proportional to the Poisson bracket {p, q}.
 851            - Higher orders include correction terms involving higher mixed derivatives.
 852    
 853        Returns
 854        -------
 855        sympy.Expr
 856            Symbolic expression for the asymptotic expansion of the commutator symbol 
 857            σ([A,B]) = σ(A∘B − B∘A).
 858    
 859        """
 860        assert self.dim == other.dim, "Operator dimensions must match"
 861        p, q = self.symbol, other.symbol
 862    
 863        pq = self.compose_asymptotic(other, order=order, mode=mode, sign_convention=sign_convention)
 864        qp = other.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
 865        
 866        comm_symbol = simplify(pq-qp)
 867
 868        return comm_symbol
 869
 870    def right_inverse_asymptotic(self, order=1):
 871        """
 872        Construct a formal right inverse R of the pseudo-differential operator P such that 
 873        the composition P ∘ R equals the identity plus a smoothing operator of order -order.
 874    
 875        This method computes an asymptotic expansion for the right inverse using recursive 
 876        corrections based on derivatives of the symbol p(x, ξ) and lower-order terms of R.
 877    
 878        Parameters
 879        ----------
 880        order : int
 881            Number of terms to include in the asymptotic expansion. Higher values improve 
 882            approximation at the cost of complexity and computational effort.
 883    
 884        Returns
 885        -------
 886        sympy.Expr
 887            The symbolic expression representing the formal right inverse R(x, ξ), which satisfies:
 888            P ∘ R = Id + O(⟨ξ⟩^{-order}), where ⟨ξ⟩ = (1 + |ξ|²)^{1/2}.
 889    
 890        Notes
 891        -----
 892        - In 1D: The recursion involves spatial derivatives of R and derivatives of p with respect to ξ.
 893        - In 2D: The multi-index generalization is used with mixed derivatives in ξ and η.
 894        - The construction relies on the non-vanishing of the principal symbol p to ensure invertibility.
 895        - Each term in the expansion corresponds to higher-order corrections involving commutators 
 896          between the operator P and the current approximation of R.
 897        """
 898        p = self.symbol
 899        if self.dim == 1:
 900            x = self.vars_x[0]
 901            xi = symbols('xi', real=True)
 902            r = 1 / p.subs(xi, xi)  # r0
 903            R = r
 904            for n in range(1, order + 1):
 905                term = 0
 906                for k in range(1, n + 1):
 907                    coeff = (1j)**(-k) / factorial(k)
 908                    inner = diff(p, xi, k) * diff(R, x, k)
 909                    term += coeff * inner
 910                R = R - r * term
 911        elif self.dim == 2:
 912            x, y = self.vars_x
 913            xi, eta = symbols('xi eta', real=True)
 914            r = 1 / p.subs({xi: xi, eta: eta})
 915            R = r
 916            for n in range(1, order + 1):
 917                term = 0
 918                for k1 in range(n + 1):
 919                    for k2 in range(n + 1 - k1):
 920                        if k1 + k2 == 0: continue
 921                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 922                        dp = diff(p, xi, k1, eta, k2)
 923                        dR = diff(R, x, k1, y, k2)
 924                        term += coeff * dp * dR
 925                R = R - r * term
 926        return R
 927
 928    def left_inverse_asymptotic(self, order=1):
 929        """
 930        Construct a formal left inverse L such that the composition L ∘ P equals the identity 
 931        operator up to terms of order ξ^{-order}. This expansion is performed asymptotically 
 932        at infinity in the frequency variable(s).
 933    
 934        The left inverse is built iteratively using symbolic differentiation and the 
 935        method of asymptotic expansions for pseudo-differential operators. It ensures that:
 936        
 937            L(P(x,ξ),x,D) ∘ P(x,D) = Id + smoothing operator of order -order
 938    
 939        Parameters
 940        ----------
 941        order : int, optional
 942            Maximum number of terms in the asymptotic expansion (default is 1). Higher values 
 943            yield more accurate inverses at the cost of increased computational complexity.
 944    
 945        Returns
 946        -------
 947        sympy.Expr
 948            Symbolic expression representing the principal symbol of the formal left inverse 
 949            operator L(x,ξ). This expression depends on spatial variables and frequencies, 
 950            and includes correction terms up to the specified order.
 951    
 952        Notes
 953        -----
 954        - In 1D: Uses recursive application of the Leibniz formula for symbols.
 955        - In 2D: Generalizes to multi-indices for mixed derivatives in (x,y) and (ξ,η).
 956        - Each term involves combinations of derivatives of the original symbol p(x,ξ) and 
 957          previously computed terms of the inverse.
 958        - Coefficients include powers of 1j (i) and factorial normalization for derivative terms.
 959        """
 960        p = self.symbol
 961        if self.dim == 1:
 962            x = self.vars_x[0]
 963            xi = symbols('xi', real=True)
 964            l = 1 / p.subs(xi, xi)
 965            L = l
 966            for n in range(1, order + 1):
 967                term = 0
 968                for k in range(1, n + 1):
 969                    coeff = (1j)**(-k) / factorial(k)
 970                    inner = diff(L, xi, k) * diff(p, x, k)
 971                    term += coeff * inner
 972                L = L - term * l
 973        elif self.dim == 2:
 974            x, y = self.vars_x
 975            xi, eta = symbols('xi eta', real=True)
 976            l = 1 / p.subs({xi: xi, eta: eta})
 977            L = l
 978            for n in range(1, order + 1):
 979                term = 0
 980                for k1 in range(n + 1):
 981                    for k2 in range(n + 1 - k1):
 982                        if k1 + k2 == 0: continue
 983                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 984                        dp = diff(p, x, k1, y, k2)
 985                        dL = diff(L, xi, k1, eta, k2)
 986                        term += coeff * dL * dp
 987                L = L - term * l
 988        return L
 989
 990    def formal_adjoint(self):
 991        """
 992        Compute the formal adjoint symbol P* of the pseudo-differential operator.
 993
 994        The adjoint is defined such that for any test functions u and v,
 995        ⟨P u, v⟩ = ⟨u, P* v⟩ holds in the distributional sense. This is obtained by 
 996        taking the complex conjugate of the symbol and expanding it asymptotically 
 997        at infinity to ensure proper behavior under integration by parts.
 998
 999        Returns
1000        -------
1001        sympy.Expr
1002            The adjoint symbol P*(x, ξ) in 1D or P*(x, y, ξ, η) in 2D.
1003        
1004        Notes:
1005        - In 1D, the expansion is performed in powers of 1/|ξ|.
1006        - In 2D, the expansion is radial in |ξ| = sqrt(ξ² + η²).
1007        - This method ensures symbolic simplifications for readability and efficiency.
1008        """
1009        p = self.symbol
1010        if self.dim == 1:
1011            x, = self.vars_x
1012            xi = symbols('xi', real=True)
1013            p_star = conjugate(p)
1014            p_star = simplify(series(p_star, xi, oo, n=6).removeO())
1015            return p_star
1016        elif self.dim == 2:
1017            x, y = self.vars_x
1018            xi, eta = symbols('xi eta', real=True)
1019            p_star = conjugate(p)
1020            p_star = simplify(series(p_star, sqrt(xi**2 + eta**2), oo, n=6).removeO())
1021            return p_star
1022
1023    def exponential_symbol(self, t=1.0, order=1, mode='kn', sign_convention=None):
1024        """
1025        Compute the symbol of exp(tP) using asymptotic expansion methods.
1026        
1027        This method calculates the exponential of a pseudo-differential operator 
1028        using either a direct power series expansion or a Magnus expansion, 
1029        depending on the structure of the symbol. The result is valid up to 
1030        the specified asymptotic order.
1031        
1032        Parameters
1033        ----------
1034        t : float or sympy.Symbol, default=1.0
1035            Time or evolution parameter. Common uses:
1036            - t = -i*τ for Schrödinger evolution: exp(-iτH)
1037            - t = τ for heat/diffusion: exp(τΔ)
1038            - t for general propagators
1039        order : int, default=3
1040            Maximum order of the asymptotic expansion. Higher orders include 
1041            more composition terms, improving accuracy for small t or when 
1042            non-commutativity effects are significant.
1043        
1044        Returns
1045        -------
1046        sympy.Expr
1047            Symbolic expression for the exponential operator symbol, computed 
1048            as an asymptotic series up to the specified order.
1049        
1050        Notes
1051        -----
1052        - For commutative symbols (e.g., pure multiplication operators), the 
1053          exponential is exact: exp(tP) = exp(t*p(x,ξ)).
1054        
1055        - For general non-commutative operators, the method uses the BCH-type 
1056          expansion via iterated composition:
1057          exp(tP) ~ I + tP + (t²/2!)P∘P + (t³/3!)P∘P∘P + ...
1058          
1059        - Each power P^n is computed via compose_asymptotic, which accounts 
1060          for the non-commutativity through derivative terms.
1061        
1062        - The expansion is valid for |t| small enough or when the symbol has 
1063          appropriate decay/growth properties.
1064        
1065        - In quantum mechanics (Schrödinger): U(t) = exp(-itH/ℏ) represents 
1066          the time evolution operator.
1067        
1068        - In parabolic PDEs (heat equation): exp(tΔ) is the heat kernel.
1069
1070        """
1071        if self.dim == 1:
1072            x = self.vars_x[0]
1073            xi = symbols('xi', real=True)
1074            
1075            # Initialize with identity
1076            result = 1
1077            
1078            # First order term: tP
1079            current_power = self.symbol
1080            result += t * current_power
1081            
1082            # Higher order terms: (t^n/n!) P^n computed via composition
1083            for n in range(2, order + 1):
1084                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1085                # We use a temporary operator for composition
1086                temp_op = PseudoDifferentialOperator(
1087                    current_power, [x], mode='symbol'
1088                )
1089                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1090                
1091                # Add term (t^n/n!) * P^n
1092                coeff = t**n / factorial(n)
1093                result += coeff * current_power
1094            
1095            return simplify(result)
1096        
1097        elif self.dim == 2:
1098            x, y = self.vars_x
1099            xi, eta = symbols('xi eta', real=True)
1100            
1101            # Initialize with identity
1102            result = 1
1103            
1104            # First order term: tP
1105            current_power = self.symbol
1106            result += t * current_power
1107            
1108            # Higher order terms: (t^n/n!) P^n computed via composition
1109            for n in range(2, order + 1):
1110                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1111                temp_op = PseudoDifferentialOperator(
1112                    current_power, [x, y], mode='symbol'
1113                )
1114                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1115                
1116                # Add term (t^n/n!) * P^n
1117                coeff = t**n / factorial(n)
1118                result += coeff * current_power
1119            
1120            return simplify(result)
1121        
1122        else:
1123            raise NotImplementedError("Only 1D and 2D operators are supported")
1124        
1125    def trace_formula(self, volume_element=None, numerical=False, 
1126                      x_bounds=None, xi_bounds=None):
1127        """
1128        Compute the semiclassical trace of the pseudo-differential operator.
1129        
1130        The trace formula relates the quantum trace of an operator to a 
1131        phase-space integral of its symbol, providing a fundamental link 
1132        between classical and quantum mechanics. This implementation supports 
1133        both symbolic and numerical integration.
1134        
1135        Parameters
1136        ----------
1137        volume_element : sympy.Expr, optional
1138            Custom volume element for the phase space integration. If None, 
1139            uses the standard Liouville measure dx dξ/(2π)^d.
1140        numerical : bool, default=False
1141            If True, perform numerical integration over specified bounds.
1142            If False, attempt symbolic integration (may fail for complex symbols).
1143        x_bounds : tuple of tuples, optional
1144            Spatial integration bounds. For 1D: ((x_min, x_max),)
1145            For 2D: ((x_min, x_max), (y_min, y_max))
1146            Required if numerical=True.
1147        xi_bounds : tuple of tuples, optional
1148            Frequency integration bounds. For 1D: ((xi_min, xi_max),)
1149            For 2D: ((xi_min, xi_max), (eta_min, eta_max))
1150            Required if numerical=True.
1151        
1152        Returns
1153        -------
1154        sympy.Expr or float
1155            The trace of the operator. Returns a symbolic expression if 
1156            numerical=False, or a float if numerical=True.
1157        
1158        Notes
1159        -----
1160        - The semiclassical trace formula states:
1161          Tr(P) = (2π)^{-d} ∫∫ p(x,ξ) dx dξ
1162          where d is the spatial dimension and p(x,ξ) is the operator symbol.
1163        
1164        - For 1D: Tr(P) = (1/2π) ∫_{-∞}^{∞} ∫_{-∞}^{∞} p(x,ξ) dx dξ
1165        
1166        - For 2D: Tr(P) = (1/4π²) ∫∫∫∫ p(x,y,ξ,η) dx dy dξ dη
1167        
1168        - This formula is exact for trace-class operators and provides an 
1169          asymptotic approximation for general pseudo-differential operators.
1170        
1171        - Physical interpretation: the trace counts the "number of states" 
1172          weighted by the observable p(x,ξ).
1173        
1174        - For projection operators (χ_Ω with χ² = χ), the trace gives the 
1175          dimension of the range, related to the phase space volume of Ω.
1176        
1177        - The factor (2π)^{-d} comes from the quantum normalization of 
1178          coherent states / Weyl quantization.
1179        """
1180        from sympy import integrate, simplify, lambdify
1181        from scipy.integrate import dblquad, nquad
1182        
1183        p = self.symbol
1184        
1185        if numerical:
1186            if x_bounds is None or xi_bounds is None:
1187                raise ValueError(
1188                    "x_bounds and xi_bounds must be provided for numerical integration"
1189                )
1190        
1191        if self.dim == 1:
1192            x, = self.vars_x
1193            xi = symbols('xi', real=True)
1194            
1195            if volume_element is None:
1196                volume_element = 1 / (2 * pi)
1197            
1198            if numerical:
1199                # Numerical integration
1200                p_func = lambdify((x, xi), p, 'numpy')
1201                (x_min, x_max), = x_bounds
1202                (xi_min, xi_max), = xi_bounds
1203                
1204                def integrand(xi_val, x_val):
1205                    return p_func(x_val, xi_val)
1206                
1207                result, error = dblquad(
1208                    integrand,
1209                    x_min, x_max,
1210                    lambda x: xi_min, lambda x: xi_max
1211                )
1212                
1213                result *= float(volume_element)
1214                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1215                return result
1216            
1217            else:
1218                # Symbolic integration
1219                integrand = p * volume_element
1220                
1221                try:
1222                    # Try to integrate over xi first, then x
1223                    integral_xi = integrate(integrand, (xi, -oo, oo))
1224                    integral_x = integrate(integral_xi, (x, -oo, oo))
1225                    return simplify(integral_x)
1226                except:
1227                    print("Warning: Symbolic integration failed. Try numerical=True")
1228                    return integrate(integrand, (xi, -oo, oo), (x, -oo, oo))
1229        
1230        elif self.dim == 2:
1231            x, y = self.vars_x
1232            xi, eta = symbols('xi eta', real=True)
1233            
1234            if volume_element is None:
1235                volume_element = 1 / (4 * pi**2)
1236            
1237            if numerical:
1238                # Numerical integration in 4D
1239                p_func = lambdify((x, y, xi, eta), p, 'numpy')
1240                (x_min, x_max), (y_min, y_max) = x_bounds
1241                (xi_min, xi_max), (eta_min, eta_max) = xi_bounds
1242                
1243                def integrand(eta_val, xi_val, y_val, x_val):
1244                    return p_func(x_val, y_val, xi_val, eta_val)
1245                
1246                result, error = nquad(
1247                    integrand,
1248                    [
1249                        [eta_min, eta_max],
1250                        [xi_min, xi_max],
1251                        [y_min, y_max],
1252                        [x_min, x_max]
1253                    ]
1254                )
1255                
1256                result *= float(volume_element)
1257                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1258                return result
1259            
1260            else:
1261                # Symbolic integration
1262                integrand = p * volume_element
1263                
1264                try:
1265                    # Integrate in order: eta, xi, y, x
1266                    integral_eta = integrate(integrand, (eta, -oo, oo))
1267                    integral_xi = integrate(integral_eta, (xi, -oo, oo))
1268                    integral_y = integrate(integral_xi, (y, -oo, oo))
1269                    integral_x = integrate(integral_y, (x, -oo, oo))
1270                    return simplify(integral_x)
1271                except:
1272                    print("Warning: Symbolic integration failed. Try numerical=True")
1273                    return integrate(
1274                        integrand,
1275                        (eta, -oo, oo), (xi, -oo, oo),
1276                        (y, -oo, oo), (x, -oo, oo)
1277                    )
1278        
1279        else:
1280            raise NotImplementedError("Only 1D and 2D operators are supported")
1281
1282    def pseudospectrum_analysis(self, x_grid, lambda_real_range, lambda_imag_range, 
1283                               epsilon_levels=[1e-1, 1e-2, 1e-3, 1e-4],
1284                               resolution=100, method='spectral', L=None, N=None):
1285        """
1286        Compute and visualize the pseudospectrum of the pseudo-differential operator.
1287        
1288        The ε-pseudospectrum is defined as:
1289            Λ_ε(A) = { λ ∈ ℂ : ‖(A - λI)^{-1}‖ ≥ ε^{-1} }
1290        
1291        This method quantizes the operator symbol into a matrix representation 
1292        and samples the resolvent norm over a grid in the complex plane.
1293        
1294        Parameters
1295        ----------
1296        x_grid : ndarray
1297            Spatial discretization grid (used if method='finite_difference')
1298        lambda_real_range : tuple
1299            Real part range of complex λ: (λ_re_min, λ_re_max)
1300        lambda_imag_range : tuple
1301            Imaginary part range: (λ_im_min, λ_im_max)
1302        epsilon_levels : list of float
1303            Contour levels for ε-pseudospectrum boundaries
1304        resolution : int
1305            Number of grid points per axis in the λ-plane
1306        method : str
1307            Discretization method:
1308            - 'spectral': FFT-based spectral differentiation (periodic, high accuracy)
1309            - 'finite_difference': Standard finite differences
1310        L : float, optional
1311            Half-domain length for spectral method (default: inferred from x_grid)
1312        N : int, optional
1313            Number of grid points for spectral method (default: len(x_grid))
1314        
1315        Returns
1316        -------
1317        dict
1318            Contains:
1319            - 'lambda_grid': meshgrid of complex λ values
1320            - 'resolvent_norm': 2D array of ‖(A - λI)^{-1}‖
1321            - 'sigma_min': 2D array of σ_min(A - λI)
1322            - 'epsilon_levels': input epsilon levels
1323            - 'eigenvalues': computed eigenvalues (if available)
1324        
1325        Notes
1326        -----
1327        - For non-self-adjoint operators, the pseudospectrum can extend far from 
1328          the actual spectrum, revealing transient behavior and non-normal dynamics.
1329        - The spectral method is preferred for smooth, periodic-like symbols.
1330        - Computational cost scales as O(resolution² × N³) due to SVD at each λ.
1331        
1332        Examples
1333        --------
1334        >>> # Analyze pseudospectrum of a non-self-adjoint operator
1335        >>> x, xi = symbols('x xi', real=True)
1336        >>> symbol = xi**2 + 1j*x*xi  # non-self-adjoint
1337        >>> op = PseudoDifferentialOperator(symbol, [x], mode='symbol')
1338        >>> result = op.pseudospectrum_analysis(
1339        ...     x_grid=np.linspace(-5, 5, 128),
1340        ...     lambda_real_range=(-2, 10),
1341        ...     lambda_imag_range=(-3, 3),
1342        ...     method='spectral'
1343        ... )
1344        """
1345        from scipy.linalg import svdvals
1346        from scipy.sparse import diags
1347        
1348        if self.dim != 1:
1349            raise NotImplementedError("Pseudospectrum analysis currently supports 1D only")
1350        
1351        # --- Step 1: Quantize the operator into a matrix ---
1352        if method == 'spectral':
1353            # Spectral (FFT) discretization
1354            if L is None:
1355                L = (x_grid[-1] - x_grid[0]) / 2.0
1356            if N is None:
1357                N = len(x_grid)
1358            
1359            x_grid_spectral = np.linspace(-L, L, N, endpoint=False)
1360            dx = x_grid_spectral[1] - x_grid_spectral[0]
1361            k = np.fft.fftfreq(N, d=dx) * 2.0 * np.pi
1362            k2 = -k**2  # symbol for -d²/dx²
1363            
1364            # Build operator matrix via spectral differentiation
1365            def apply_operator(u):
1366                """Apply Op(symbol) to vector u"""
1367                u_hat = np.fft.fft(u)
1368                # Extract kinetic part from symbol (assuming symbol = f(xi) + g(x))
1369                # This is a simplified model; for general symbols, use full quantization
1370                kinetic = k2 * u_hat
1371                v = np.fft.ifft(kinetic)
1372                # Add potential/position-dependent part
1373                x_vals = x_grid_spectral
1374                potential = self.p_func(x_vals, 0.0)  # evaluate at ξ=0 for position part
1375                v += potential * u
1376                return np.real(v)
1377            
1378            # Assemble matrix
1379            H = np.zeros((N, N), dtype=complex)
1380            for j in range(N):
1381                e = np.zeros(N)
1382                e[j] = 1.0
1383                H[:, j] = apply_operator(e)
1384            
1385            print(f"Operator quantized via spectral method: {N}×{N} matrix")
1386        
1387        elif method == 'finite_difference':
1388            # Finite difference discretization
1389            N = len(x_grid)
1390            dx = x_grid[1] - x_grid[0]
1391            
1392            # Build -d²/dx² using centered differences
1393            diag_main = -2.0 / dx**2 * np.ones(N)
1394            diag_off = 1.0 / dx**2 * np.ones(N - 1)
1395            D2 = diags([diag_off, diag_main, diag_off], [-1, 0, 1], shape=(N, N)).toarray()
1396            
1397            # Add position-dependent part from symbol
1398            x_vals = x_grid
1399            potential = np.diag(self.p_func(x_vals, 0.0))
1400            
1401            H = -D2 + potential
1402            print(f"Operator quantized via finite differences: {N}×{N} matrix")
1403        
1404        else:
1405            raise ValueError("method must be 'spectral' or 'finite_difference'")
1406        
1407        # --- Step 2: Sample resolvent norm over λ-plane ---
1408        lambda_re = np.linspace(*lambda_real_range, resolution)
1409        lambda_im = np.linspace(*lambda_imag_range, resolution)
1410        Lambda_re, Lambda_im = np.meshgrid(lambda_re, lambda_im)
1411        Lambda = Lambda_re + 1j * Lambda_im
1412        
1413        resolvent_norm = np.zeros_like(Lambda, dtype=float)
1414        sigma_min_grid = np.zeros_like(Lambda, dtype=float)
1415        
1416        I = np.eye(N)
1417        
1418        print(f"Computing pseudospectrum over {resolution}×{resolution} grid...")
1419        for i in range(resolution):
1420            for j in range(resolution):
1421                lam = Lambda[i, j]
1422                A = H - lam * I
1423                
1424                try:
1425                    # Compute smallest singular value
1426                    s = svdvals(A)
1427                    s_min = s[-1]
1428                    sigma_min_grid[i, j] = s_min
1429                    resolvent_norm[i, j] = 1.0 / (s_min + 1e-16)  # regularization
1430                except Exception:
1431                    resolvent_norm[i, j] = np.nan
1432                    sigma_min_grid[i, j] = np.nan
1433        
1434        # --- Step 3: Compute eigenvalues ---
1435        try:
1436            eigenvalues = np.linalg.eigvals(H)
1437        except:
1438            eigenvalues = None
1439        
1440        # --- Step 4: Visualization ---
1441        plt.figure(figsize=(14, 6))
1442        
1443        # Left panel: log10(resolvent norm)
1444        plt.subplot(1, 2, 1)
1445        levels_log = np.log10(1.0 / np.array(epsilon_levels))
1446        cs = plt.contour(Lambda_re, Lambda_im, np.log10(resolvent_norm + 1e-16), 
1447                         levels=levels_log, colors='blue', linewidths=1.5)
1448        plt.clabel(cs, inline=True, fmt='ε=10^%d')
1449        
1450        if eigenvalues is not None:
1451            plt.plot(eigenvalues.real, eigenvalues.imag, 'r*', markersize=8, label='Eigenvalues')
1452        
1453        plt.xlabel('Re(λ)')
1454        plt.ylabel('Im(λ)')
1455        plt.title('ε-Pseudospectrum: log₁₀(‖(A - λI)⁻¹‖)')
1456        plt.grid(alpha=0.3)
1457        plt.legend()
1458        plt.axis('equal')
1459        
1460        # Right panel: σ_min contours
1461        plt.subplot(1, 2, 2)
1462        cs2 = plt.contourf(Lambda_re, Lambda_im, sigma_min_grid, 
1463                           levels=50, cmap='viridis')
1464        plt.colorbar(cs2, label='σ_min(A - λI)')
1465        
1466        if eigenvalues is not None:
1467            plt.plot(eigenvalues.real, eigenvalues.imag, 'r*', markersize=8)
1468        
1469        for eps in epsilon_levels:
1470            plt.contour(Lambda_re, Lambda_im, sigma_min_grid, 
1471                       levels=[eps], colors='red', linewidths=1.5, alpha=0.7)
1472        
1473        plt.xlabel('Re(λ)')
1474        plt.ylabel('Im(λ)')
1475        plt.title('Smallest singular value σ_min(A - λI)')
1476        plt.grid(alpha=0.3)
1477        plt.axis('equal')
1478        
1479        plt.tight_layout()
1480        plt.show()
1481        
1482        return {
1483            'lambda_grid': Lambda,
1484            'resolvent_norm': resolvent_norm,
1485            'sigma_min': sigma_min_grid,
1486            'epsilon_levels': epsilon_levels,
1487            'eigenvalues': eigenvalues,
1488            'operator_matrix': H
1489        }
1490    
1491    def symplectic_flow(self):
1492        """
1493        Compute the Hamiltonian vector field associated with the principal symbol.
1494
1495        This method derives the canonical equations of motion for the phase space variables 
1496        (x, ξ) in 1D or (x, y, ξ, η) in 2D, based on the Hamiltonian formalism. These describe 
1497        how position and frequency variables evolve under the flow generated by the symbol.
1498
1499        Returns
1500        -------
1501        dict
1502            A dictionary containing the components of the Hamiltonian vector field:
1503            - In 1D: keys are 'dx/dt' and 'dxi/dt', corresponding to dx/dt = ∂p/∂ξ and dξ/dt = -∂p/∂x.
1504            - In 2D: keys are 'dx/dt', 'dy/dt', 'dxi/dt', and 'deta/dt', with similar definitions:
1505              dx/dt = ∂p/∂ξ, dy/dt = ∂p/∂η, dξ/dt = -∂p/∂x, dη/dt = -∂p/∂y.
1506
1507        Notes
1508        -----
1509        - The Hamiltonian here is the principal symbol p(x, ξ) itself.
1510        - This flow preserves the symplectic structure of phase space.
1511        """
1512        if self.dim == 1:
1513            x,  = self.vars_x
1514            xi = symbols('xi', real=True)
1515            return {
1516                'dx/dt': diff(self.symbol, xi),
1517                'dxi/dt': -diff(self.symbol, x)
1518            }
1519        elif self.dim == 2:
1520            x, y = self.vars_x
1521            xi, eta = symbols('xi eta', real=True)
1522            return {
1523                'dx/dt': diff(self.symbol, xi),
1524                'dy/dt': diff(self.symbol, eta),
1525                'dxi/dt': -diff(self.symbol, x),
1526                'deta/dt': -diff(self.symbol, y)
1527            }
1528
1529    def is_elliptic_numerically(self, x_grid, xi_grid, threshold=1e-8):
1530        """
1531        Check if the pseudo-differential symbol p(x, ξ) is elliptic over a given grid.
1532    
1533        A symbol is considered elliptic if its magnitude |p(x, ξ)| remains bounded away from zero 
1534        across all points in the spatial-frequency domain. This method evaluates the symbol on a 
1535        grid of spatial and frequency coordinates and checks whether its minimum absolute value 
1536        exceeds a specified threshold.
1537    
1538        Resampling is applied to large grids to prevent excessive memory usage, particularly in 2D.
1539    
1540        Parameters
1541        ----------
1542        x_grid : ndarray
1543            Spatial grid: either a 1D array (x) or a tuple of two 1D arrays (x, y).
1544        xi_grid : ndarray
1545            Frequency grid: either a 1D array (ξ) or a tuple of two 1D arrays (ξ, η).
1546        threshold : float, optional
1547            Minimum acceptable value for |p(x, ξ)|. If the smallest evaluated symbol value falls below this,
1548            the symbol is not considered elliptic.
1549    
1550        Returns
1551        -------
1552        bool
1553            True if the symbol is elliptic on the resampled grid, False otherwise.
1554        """
1555        RESAMPLE_SIZE = 32  # Reduced size to prevent memory explosion
1556        
1557        if self.dim == 1:
1558            x_vals = x_grid
1559            xi_vals = xi_grid
1560            # Resampling if necessary
1561            if len(x_vals) > RESAMPLE_SIZE:
1562                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1563            if len(xi_vals) > RESAMPLE_SIZE:
1564                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1565        
1566            X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1567            symbol_vals = self.p_func(X, XI)
1568        
1569        elif self.dim == 2:
1570            x_vals, y_vals = x_grid
1571            xi_vals, eta_vals = xi_grid
1572        
1573            # Spatial resampling
1574            if len(x_vals) > RESAMPLE_SIZE:
1575                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1576            if len(y_vals) > RESAMPLE_SIZE:
1577                y_vals = np.linspace(y_vals.min(), y_vals.max(), RESAMPLE_SIZE)
1578        
1579            # Frequency resampling
1580            if len(xi_vals) > RESAMPLE_SIZE:
1581                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1582            if len(eta_vals) > RESAMPLE_SIZE:
1583                eta_vals = np.linspace(eta_vals.min(), eta_vals.max(), RESAMPLE_SIZE)
1584        
1585            X, Y, XI, ETA = np.meshgrid(x_vals, y_vals, xi_vals, eta_vals, indexing='ij')
1586            symbol_vals = self.p_func(X, Y, XI, ETA)
1587        
1588        min_abs_val = np.min(np.abs(symbol_vals))
1589        return min_abs_val > threshold
1590
1591
1592    def is_self_adjoint(self, tol=1e-10):
1593        """
1594        Check whether the pseudo-differential operator is formally self-adjoint (Hermitian).
1595
1596        A self-adjoint operator satisfies P = P*, where P* is the formal adjoint of P.
1597        This property is essential for ensuring real-valued eigenvalues and stable evolution 
1598        in quantum mechanics and symmetric wave propagation.
1599
1600        Parameters
1601        ----------
1602        tol : float
1603            Tolerance for symbolic comparison between P and P*. Small numerical differences 
1604            below this threshold are considered equal.
1605
1606        Returns
1607        -------
1608        bool
1609            True if the symbol p(x, ξ) equals its formal adjoint p*(x, ξ) within the given tolerance,
1610            indicating that the operator is self-adjoint.
1611
1612        Notes:
1613        - The formal adjoint is computed via conjugation and asymptotic expansion at infinity in ξ.
1614        - Symbolic simplification is used to verify equality, ensuring robustness against superficial 
1615          expression differences.
1616        """
1617        p = self.symbol
1618        p_star = self.formal_adjoint()
1619        return simplify(p - p_star).equals(0)
1620
1621    def visualize_fiber(self, x_grid, xi_grid, x0=0.0, y0=0.0):
1622        """
1623        Plot the cotangent fiber structure at a fixed spatial point (x₀[, y₀]).
1624    
1625        This visualization shows how the symbol p(x, ξ) behaves on the cotangent fiber 
1626        above a fixed spatial point. In microlocal analysis, this provides insight into 
1627        the frequency content of the operator at that location.
1628    
1629        Parameters
1630        ----------
1631        x_grid : ndarray
1632            Spatial grid values (1D) for evaluation in 1D case.
1633        xi_grid : ndarray
1634            Frequency grid values (1D) for evaluation in both 1D and 2D cases.
1635        x0 : float, optional
1636            Fixed x-coordinate of the base point in space (1D or 2D).
1637        y0 : float, optional
1638            Fixed y-coordinate of the base point in space (2D only).
1639    
1640        Notes
1641        -----
1642        - In 1D: Displays |p(x, ξ)| over the (x, ξ) phase plane near the fixed point.
1643        - In 2D: Fixes (x₀, y₀) and evaluates p(x₀, y₀, ξ, η), showing the fiber over that point.
1644        - The color map represents the magnitude of the symbol, highlighting regions where it vanishes or becomes singular.
1645    
1646        Raises
1647        ------
1648        NotImplementedError
1649            If called in 2D with missing or improperly formatted grids.
1650        """
1651        if self.dim == 1:
1652            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1653            symbol_vals = self.p_func(X, XI)
1654            plt.contourf(X, XI, np.abs(symbol_vals), levels=50, cmap='viridis')
1655            plt.colorbar(label='|Symbol|')
1656            plt.xlabel('x (position)')
1657            plt.ylabel('ξ (frequency)')
1658            plt.title('Cotangent Fiber Structure')
1659            plt.show()
1660        elif self.dim == 2:
1661            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, xi_grid)
1662            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1663            plt.contourf(xi_grid, xi_grid, np.abs(symbol_vals), levels=50, cmap='viridis')
1664            plt.colorbar(label='|Symbol|')
1665            plt.xlabel('ξ')
1666            plt.ylabel('η')
1667            plt.title(f'Cotangent Fiber at x={x0}, y={y0}')
1668            plt.show()
1669
1670    def visualize_symbol_amplitude(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
1671        """
1672        Display the modulus |p(x, ξ)| or |p(x, y, ξ₀, η₀)| as a color map.
1673    
1674        This method visualizes the amplitude of the pseudodifferential operator's symbol 
1675        in either 1D or 2D spatial configuration. In 2D, the frequency variables are fixed 
1676        to specified values (ξ₀, η₀) for visualization purposes.
1677    
1678        Parameters
1679        ----------
1680        x_grid, y_grid : ndarray
1681            Spatial grids over which to evaluate the symbol. y_grid is optional and used only in 2D.
1682        xi_grid, eta_grid : ndarray
1683            Frequency grids. In 2D, these define the domain over which the symbol is evaluated,
1684            but the visualization fixes ξ = ξ₀ and η = η₀.
1685        xi0, eta0 : float, optional
1686            Fixed frequency values for slicing in 2D visualization. Defaults to zero.
1687    
1688        Notes
1689        -----
1690        - In 1D: Visualizes |p(x, ξ)| over the (x, ξ) grid.
1691        - In 2D: Visualizes |p(x, y, ξ₀, η₀)| at fixed frequencies ξ₀ and η₀.
1692        - The color intensity represents the magnitude of the symbol, highlighting regions where the symbol is large or small.
1693        """
1694        if self.dim == 1:
1695            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1696            symbol_vals = self.p_func(X, XI) 
1697            plt.pcolormesh(X, XI, np.abs(symbol_vals), shading='auto')
1698            plt.colorbar(label='|Symbol|')
1699            plt.xlabel('x')
1700            plt.ylabel('ξ')
1701            plt.title('Symbol Amplitude |p(x, ξ)|')
1702            plt.show()
1703        elif self.dim == 2:
1704            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
1705            XI = np.full_like(X, xi0)
1706            ETA = np.full_like(Y, eta0)
1707            symbol_vals = self.p_func(X, Y, XI, ETA)
1708            plt.pcolormesh(X, Y, np.abs(symbol_vals), shading='auto')
1709            plt.colorbar(label='|Symbol|')
1710            plt.xlabel('x')
1711            plt.ylabel('y')
1712            plt.title(f'Symbol Amplitude at ξ={xi0}, η={eta0}')
1713            plt.show()
1714
1715    def visualize_phase(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
1716        """
1717        Plot the phase (argument) of the pseudodifferential operator's symbol p(x, ξ) or p(x, y, ξ, η).
1718
1719        This visualization helps in understanding the oscillatory behavior and regularity 
1720        properties of the operator in phase space. The phase is displayed modulo 2π using 
1721        a cyclic colormap ('twilight') to emphasize its periodic nature.
1722
1723        Parameters
1724        ----------
1725        x_grid : ndarray
1726            1D array of spatial coordinates (x).
1727        xi_grid : ndarray
1728            1D array of frequency coordinates (ξ).
1729        y_grid : ndarray, optional
1730            2D spatial grid for y-coordinate (in 2D problems). Default is None.
1731        eta_grid : ndarray, optional
1732            2D frequency grid for η (in 2D problems). Not used directly but kept for API consistency.
1733        xi0 : float, optional
1734            Fixed value of ξ for slicing in 2D visualization. Default is 0.0.
1735        eta0 : float, optional
1736            Fixed value of η for slicing in 2D visualization. Default is 0.0.
1737
1738        Notes:
1739        - In 1D: Displays arg(p(x, ξ)) over the (x, ξ) phase plane.
1740        - In 2D: Displays arg(p(x, y, ξ₀, η₀)) for fixed frequency values (ξ₀, η₀).
1741        - Uses plt.pcolormesh with 'twilight' colormap to represent angles from -π to π.
1742
1743        Raises:
1744        - NotImplementedError: If the spatial dimension is not 1D or 2D.
1745        """
1746        if self.dim == 1:
1747            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1748            symbol_vals = self.p_func(X, XI) 
1749            plt.pcolormesh(X, XI, np.angle(symbol_vals), shading='auto', cmap='twilight')
1750            plt.colorbar(label='arg(Symbol) [rad]')
1751            plt.xlabel('x')
1752            plt.ylabel('ξ')
1753            plt.title('Phase Portrait (arg p(x, ξ))')
1754            plt.show()
1755        elif self.dim == 2:
1756            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
1757            XI = np.full_like(X, xi0)
1758            ETA = np.full_like(Y, eta0)
1759            symbol_vals = self.p_func(X, Y, XI, ETA)
1760            plt.pcolormesh(X, Y, np.angle(symbol_vals), shading='auto', cmap='twilight')
1761            plt.colorbar(label='arg(Symbol) [rad]')
1762            plt.xlabel('x')
1763            plt.ylabel('y')
1764            plt.title(f'Phase Portrait at ξ={xi0}, η={eta0}')
1765            plt.show()
1766            
1767    def visualize_characteristic_set(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0, levels=[1e-1]):
1768        """
1769        Visualize the characteristic set of the pseudo-differential symbol, defined as the approximate zero set p(x, ξ) ≈ 0.
1770    
1771        In microlocal analysis, the characteristic set is the locus of points in phase space (x, ξ) where the symbol p(x, ξ) vanishes,
1772        playing a key role in understanding propagation of singularities.
1773    
1774        Parameters
1775        ----------
1776        x_grid : ndarray
1777            Spatial grid values (1D array) for plotting in 1D or evaluation point in 2D.
1778        xi_grid : ndarray
1779            Frequency variable grid values (1D array) used to construct the frequency domain.
1780        x0 : float, optional
1781            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific x position.
1782        y0 : float, optional
1783            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific y position.
1784    
1785        Notes
1786        -----
1787        - For 1D, this method plots the contour of |p(x, ξ)| = ε with ε = 1e-5 over the (x, ξ) plane.
1788        - For 2D, it evaluates the symbol at fixed (x₀, y₀) and plots the characteristic set in the (ξ, η) frequency plane.
1789        - This visualization helps identify directions of degeneracy or hypoellipticity of the operator.
1790    
1791        Raises
1792        ------
1793        NotImplementedError
1794            If called on a solver with dimensionality other than 1D or 2D.
1795    
1796        Displays
1797        ------
1798        A matplotlib contour plot showing either:
1799            - The characteristic curve in the (x, ξ) phase plane (1D),
1800            - The characteristic surface slice in the (ξ, η) frequency plane at (x₀, y₀) (2D).
1801        """
1802        if self.dim == 1:
1803            x_grid = np.asarray(x_grid)
1804            xi_grid = np.asarray(xi_grid)
1805            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1806            symbol_vals = self.p_func(X, XI) 
1807            plt.contour(X, XI, np.abs(symbol_vals), levels=levels, colors='red')
1808            plt.xlabel('x')
1809            plt.ylabel('ξ')
1810            plt.title('Characteristic Set (p(x, ξ) ≈ 0)')
1811            plt.grid(True)
1812            plt.show()
1813        elif self.dim == 2:
1814            if eta_grid is None:
1815                raise ValueError("eta_grid must be provided for 2D visualization.")
1816            xi_grid = np.asarray(xi_grid)
1817            eta_grid = np.asarray(eta_grid)
1818            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
1819            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1820            plt.contour(xi_grid, eta_grid, np.abs(symbol_vals), levels=levels, colors='red')
1821            plt.xlabel('ξ')
1822            plt.ylabel('η')
1823            plt.title(f'Characteristic Set at x={x0}, y={y0}')
1824            plt.grid(True)
1825            plt.show()
1826        else:
1827            raise NotImplementedError("Only 1D/2D characteristic sets supported.")
1828
1829    def visualize_characteristic_gradient(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0):
1830        """
1831        Visualize the norm of the gradient of the symbol in phase space.
1832        
1833        This method computes the magnitude of the gradient |∇p| of a pseudo-differential 
1834        symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. The resulting colormap reveals 
1835        regions where the symbol varies rapidly or remains nearly stationary, 
1836        which is particularly useful for analyzing characteristic sets.
1837        
1838        Parameters
1839        ----------
1840        x_grid : numpy.ndarray
1841            1D array of spatial coordinates for the x-direction.
1842        xi_grid : numpy.ndarray
1843            1D array of frequency coordinates (ξ).
1844        y_grid : numpy.ndarray, optional
1845            1D array of spatial coordinates for the y-direction (used in 2D mode). Default is None.
1846        eta_grid : numpy.ndarray, optional
1847            1D array of frequency coordinates (η) for the 2D case. Default is None.
1848        x0 : float, optional
1849            Fixed x-coordinate for evaluating the symbol in 2D. Default is 0.0.
1850        y0 : float, optional
1851            Fixed y-coordinate for evaluating the symbol in 2D. Default is 0.0.
1852        
1853        Returns
1854        -------
1855        None
1856            Displays a 2D colormap of |∇p| over the relevant phase-space domain.
1857        
1858        Notes
1859        -----
1860        - In 1D, the full gradient ∇p = (∂ₓp, ∂ξp) is computed over the (x, ξ) grid.
1861        - In 2D, the gradient ∇p = (∂ξp, ∂ηp) is computed at a fixed spatial point (x₀, y₀) over the (ξ, η) grid.
1862        - Numerical differentiation is performed using `np.gradient`.
1863        - High values of |∇p| indicate rapid variation of the symbol, while low values typically suggest characteristic regions.
1864        """
1865        if self.dim == 1:
1866            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1867            symbol_vals = self.p_func(X, XI)
1868            grad_x = np.gradient(symbol_vals, axis=0)
1869            grad_xi = np.gradient(symbol_vals, axis=1)
1870            grad_norm = np.sqrt(grad_x**2 + grad_xi**2)
1871            plt.pcolormesh(X, XI, grad_norm, cmap='inferno', shading='auto')
1872            plt.colorbar(label='|∇p|')
1873            plt.xlabel('x')
1874            plt.ylabel('ξ')
1875            plt.title('Gradient Norm (High Near Zeros)')
1876            plt.grid(True)
1877            plt.show()
1878        elif self.dim == 2:
1879            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
1880            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1881            grad_xi = np.gradient(symbol_vals, axis=0)
1882            grad_eta = np.gradient(symbol_vals, axis=1)
1883            grad_norm = np.sqrt(np.abs(grad_xi)**2 + np.abs(grad_eta)**2)
1884            plt.pcolormesh(xi_grid, eta_grid, grad_norm, cmap='inferno', shading='auto')
1885            plt.colorbar(label='|∇p|')
1886            plt.xlabel('ξ')
1887            plt.ylabel('η')
1888            plt.title(f'Gradient Norm at x={x0}, y={y0}')
1889            plt.grid(True)
1890            plt.show()
1891
1892    def plot_hamiltonian_flow(self, x0=0.0, xi0=5.0, y0=0.0, eta0=0.0, tmax=1.0, n_steps=100, show_field=True):
1893        """
1894        Integrate and plot the Hamiltonian trajectories of the symbol in phase space.
1895
1896        This method numerically integrates the Hamiltonian vector field derived from 
1897        the operator's symbol to visualize how singularities propagate under the flow. 
1898        It supports both 1D and 2D problems.
1899
1900        Parameters
1901        ----------
1902        x0, xi0 : float
1903            Initial position and frequency (momentum) in 1D.
1904        y0, eta0 : float, optional
1905            Initial position and frequency in 2D; defaults to zero.
1906        tmax : float
1907            Final integration time for the ODE solver.
1908        n_steps : int
1909            Number of time steps used in the integration.
1910
1911        Notes
1912        -----
1913        - The Hamiltonian vector field is obtained from the symplectic flow of the symbol.
1914        - If the field is complex-valued, only its real part is used for integration.
1915        - In 1D, the trajectory is plotted in (x, ξ) phase space.
1916        - In 2D, the spatial trajectory (x(t), y(t)) is shown along with instantaneous 
1917          momentum vectors (ξ(t), η(t)) using a quiver plot.
1918
1919        Raises
1920        ------
1921        NotImplementedError
1922            If the spatial dimension is not 1D or 2D.
1923
1924        Displays
1925        --------
1926        matplotlib plot
1927            Phase space trajectory(ies) showing the evolution of position and momentum 
1928            under the Hamiltonian dynamics.
1929        """
1930        def make_real(expr):
1931            from sympy import re, simplify
1932            expr = expr.doit(deep=True)
1933            return simplify(re(expr))
1934    
1935        H = self.symplectic_flow()
1936    
1937        if any(im(H[k]) != 0 for k in H):
1938            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
1939    
1940        if self.dim == 1:
1941            x, = self.vars_x
1942            xi = symbols('xi', real=True)
1943    
1944            dxdt_expr = make_real(H['dx/dt'])
1945            dxidt_expr = make_real(H['dxi/dt'])
1946    
1947            dxdt = lambdify((x, xi), dxdt_expr, 'numpy')
1948            dxidt = lambdify((x, xi), dxidt_expr, 'numpy')
1949    
1950            def hamilton(t, Y):
1951                x, xi = Y
1952                return [dxdt(x, xi), dxidt(x, xi)]
1953    
1954            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0], t_eval=np.linspace(0, tmax, n_steps))
1955
1956            if sol.status != 0:
1957                print(f"⚠️ Integration warning: {sol.message}")
1958            
1959            n_points = sol.y.shape[1]
1960            if n_points < n_steps:
1961                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
1962                n_steps = n_points
1963
1964            x_vals, xi_vals = sol.y
1965    
1966            plt.plot(x_vals, xi_vals)
1967            plt.xlabel("x")
1968            plt.ylabel("ξ")
1969            plt.title("Hamiltonian Flow in Phase Space (1D)")
1970            plt.grid(True)
1971            plt.show()
1972    
1973        elif self.dim == 2:
1974            x, y = self.vars_x
1975            xi, eta = symbols('xi eta', real=True)
1976    
1977            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
1978            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
1979            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
1980            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
1981    
1982            def hamilton(t, Y):
1983                x, y, xi, eta = Y
1984                return [
1985                    dxdt(x, y, xi, eta),
1986                    dydt(x, y, xi, eta),
1987                    dxidt(x, y, xi, eta),
1988                    detadt(x, y, xi, eta)
1989                ]
1990    
1991            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0], t_eval=np.linspace(0, tmax, n_steps))
1992
1993            if sol.status != 0:
1994                print(f"⚠️ Integration warning: {sol.message}")
1995            
1996            n_points = sol.y.shape[1]
1997            if n_points < n_steps:
1998                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
1999                n_steps = n_points
2000
2001            x_vals, y_vals, xi_vals, eta_vals = sol.y
2002    
2003            plt.plot(x_vals, y_vals, label='Position')
2004            plt.quiver(x_vals, y_vals, xi_vals, eta_vals, scale=20, width=0.003, alpha=0.5, color='r')
2005            
2006            # Vector field of the flow (optional)
2007            if show_field:
2008                X, Y = np.meshgrid(np.linspace(min(x_vals), max(x_vals), 20),
2009                                   np.linspace(min(y_vals), max(y_vals), 20))
2010                XI, ETA = xi0 * np.ones_like(X), eta0 * np.ones_like(Y)
2011                U = dxdt(X, Y, XI, ETA)
2012                V = dydt(X, Y, XI, ETA)
2013                plt.quiver(X, Y, U, V, color='gray', alpha=0.2, scale=30, width=0.002)
2014
2015            plt.xlabel("x")
2016            plt.ylabel("y")
2017            plt.title("Hamiltonian Flow in Phase Space (2D)")
2018            plt.legend()
2019            plt.grid(True)
2020            plt.axis('equal')
2021            plt.show()
2022
2023    def plot_symplectic_vector_field(self, xlim=(-2, 2), klim=(-5, 5), density=30):
2024        """
2025        Visualize the symplectic vector field (Hamiltonian vector field) associated with the operator's symbol.
2026
2027        The plotted vector field corresponds to (∂_ξ p, -∂_x p), where p(x, ξ) is the principal symbol 
2028        of the pseudo-differential operator. This field governs the bicharacteristic flow in phase space.
2029
2030        Parameters
2031        ----------
2032        xlim : tuple of float
2033            Range for spatial variable x, as (x_min, x_max).
2034        klim : tuple of float
2035            Range for frequency variable ξ, as (ξ_min, ξ_max).
2036        density : int
2037            Number of grid points per axis for the visualization grid.
2038
2039        Raises
2040        ------
2041        NotImplementedError
2042            If called on a 2D operator (currently only 1D implementation available).
2043
2044        Notes
2045        -----
2046        - Only supports one-dimensional operators.
2047        - Uses symbolic differentiation to compute ∂_ξ p and ∂_x p.
2048        - Numerical evaluation is done via lambdify with NumPy backend.
2049        - Visualization uses matplotlib quiver plot to show vector directions.
2050        """
2051        x_vals = np.linspace(*xlim, density)
2052        xi_vals = np.linspace(*klim, density)
2053        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2054
2055        if self.dim != 1:
2056            raise NotImplementedError("Only 1D version implemented.")
2057
2058        x, = self.vars_x
2059        xi = symbols('xi', real=True)
2060        H = self.symplectic_flow()
2061        dxdt = lambdify((x, xi), simplify(H['dx/dt']), 'numpy')
2062        dxidt = lambdify((x, xi), simplify(H['dxi/dt']), 'numpy')
2063
2064        U = dxdt(X, XI)
2065        V = dxidt(X, XI)
2066
2067        plt.quiver(X, XI, U, V, scale=10, width=0.005)
2068        plt.xlabel('x')
2069        plt.ylabel(r'$\xi$')
2070        plt.title("Symplectic Vector Field (1D)")
2071        plt.grid(True)
2072        plt.show()
2073
2074    def visualize_micro_support(self, xlim=(-2, 2), klim=(-10, 10), threshold=1e-3, density=300):
2075        """
2076        Visualize the micro-support of the operator by plotting the inverse of the symbol magnitude 1 / |p(x, ξ)|.
2077    
2078        The micro-support provides insight into the singularities of a pseudo-differential operator 
2079        in phase space (x, ξ). Regions where |p(x, ξ)| is small correspond to large values in 1/|p(x, ξ)|,
2080        highlighting areas of significant operator influence or singularity.
2081    
2082        Parameters
2083        ----------
2084        xlim : tuple
2085            Spatial domain limits (x_min, x_max).
2086        klim : tuple
2087            Frequency domain limits (ξ_min, ξ_max).
2088        threshold : float
2089            Threshold below which |p(x, ξ)| is considered effectively zero; used for numerical stability.
2090        density : int
2091            Number of grid points along each axis for visualization resolution.
2092    
2093        Raises
2094        ------
2095        NotImplementedError
2096            If called on a solver with dimension greater than 1 (only 1D visualization is supported).
2097    
2098        Notes
2099        -----
2100        - This method evaluates the symbol p(x, ξ) over a grid and plots its reciprocal to emphasize 
2101          regions where the symbol is near zero.
2102        - A small constant (1e-10) is added to the denominator to avoid division by zero.
2103        - The resulting plot helps identify characteristic sets.
2104        """
2105        if self.dim != 1:
2106            raise NotImplementedError("Only 1D micro-support visualization implemented.")
2107
2108        x_vals = np.linspace(*xlim, density)
2109        xi_vals = np.linspace(*klim, density)
2110        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2111        Z = np.abs(self.p_func(X, XI))
2112
2113        plt.contourf(X, XI, 1 / (Z + 1e-10), levels=100, cmap='inferno')
2114        plt.colorbar(label=r'$1/|p(x,\xi)|$')
2115        plt.xlabel('x')
2116        plt.ylabel(r'$\xi$')
2117        plt.title("Micro-Support Estimate (1/|Symbol|)")
2118        plt.show()
2119
2120    def group_velocity_field(self, xlim=(-2, 2), klim=(-10, 10), density=30):
2121        """
2122        Plot the group velocity field ∇_ξ p(x, ξ) for 1D pseudo-differential operators.
2123
2124        The group velocity represents the speed at which waves of different frequencies propagate 
2125        in a dispersive medium. It is defined as the gradient of the symbol p(x, ξ) with respect 
2126        to the frequency variable ξ.
2127
2128        Parameters
2129        ----------
2130        xlim : tuple of float
2131            Spatial domain limits (x-axis).
2132        klim : tuple of float
2133            Frequency domain limits (ξ-axis).
2134        density : int
2135            Number of grid points per axis used for visualization.
2136
2137        Raises
2138        ------
2139        NotImplementedError
2140            If called on a 2D operator, since this visualization is only implemented for 1D.
2141
2142        Notes
2143        -----
2144        - This method visualizes the vector field (∂p/∂ξ) in phase space.
2145        - Used for analyzing wave propagation properties and dispersion relations.
2146        - Requires symbolic expression self.expr depending on x and ξ.
2147        """
2148        if self.dim != 1:
2149            raise NotImplementedError("Only 1D group velocity visualization implemented.")
2150
2151        x, = self.vars_x
2152        xi = symbols('xi', real=True)
2153        dp_dxi = diff(self.symbol, xi)
2154        grad_func = lambdify((x, xi), dp_dxi, 'numpy')
2155
2156        x_vals = np.linspace(*xlim, density)
2157        xi_vals = np.linspace(*klim, density)
2158        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2159        V = grad_func(X, XI)
2160
2161        plt.quiver(X, XI, np.ones_like(V), V, scale=10, width=0.004)
2162        plt.xlabel('x')
2163        plt.ylabel(r'$\xi$')
2164        plt.title("Group Velocity Field (1D)")
2165        plt.grid(True)
2166        plt.show()
2167
2168    def animate_singularity(self, xi0=5.0, eta0=0.0, x0=0.0, y0=0.0,
2169                            tmax=4.0, n_frames=100, projection=None):
2170        """
2171        Animate the propagation of a singularity under the Hamiltonian flow.
2172
2173        This method visualizes how a singularity (x₀, y₀, ξ₀, η₀) evolves in phase space 
2174        according to the Hamiltonian dynamics induced by the principal symbol of the operator.
2175        The animation integrates the Hamiltonian equations of motion and supports various projections:
2176        position (x-y), frequency (ξ-η), or mixed phase space coordinates.
2177
2178        Parameters
2179        ----------
2180        xi0, eta0 : float
2181            Initial frequency components (ξ₀, η₀).
2182        x0, y0 : float
2183            Initial spatial coordinates (x₀, y₀).
2184        tmax : float
2185            Total time of integration (final animation time).
2186        n_frames : int
2187            Number of frames in the resulting animation.
2188        projection : str or None
2189            Type of projection to display:
2190                - 'position' : x vs y (or x alone in 1D)
2191                - 'frequency': ξ vs η (or ξ alone in 1D)
2192                - 'phase'    : mixed coordinates like x vs ξ or x vs η
2193                If None, defaults to 'phase' in 1D and 'position' in 2D.
2194
2195        Returns
2196        -------
2197        matplotlib.animation.FuncAnimation
2198            Animation object that can be displayed interactively in Jupyter notebooks or saved as a video.
2199
2200        Notes
2201        -----
2202        - In 1D, only one spatial and one frequency variable are used.
2203        - Complex-valued Hamiltonian fields are truncated to their real parts for integration.
2204        - Trajectories are shown with both instantaneous position (dot) and full path (dashed line).
2205        """
2206        rc('animation', html='jshtml')
2207    
2208        def make_real(expr):
2209            from sympy import re, simplify
2210            expr = expr.doit(deep=True)
2211            return simplify(re(expr))
2212  
2213        H = self.symplectic_flow()
2214
2215        H = {k: v.doit(deep=True) for k, v in H.items()}
2216
2217        print("H = ", H)
2218    
2219        if any(im(H[k]) != 0 for k in H):
2220            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
2221    
2222        if self.dim == 1:
2223            x, = self.vars_x
2224            xi = symbols('xi', real=True)
2225    
2226            dxdt = lambdify((x, xi), make_real(H['dx/dt']), 'numpy')
2227            dxidt = lambdify((x, xi), make_real(H['dxi/dt']), 'numpy')
2228    
2229            def hamilton(t, Y):
2230                x, xi = Y
2231                return [dxdt(x, xi), dxidt(x, xi)]
2232    
2233            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0],
2234                            t_eval=np.linspace(0, tmax, n_frames))
2235            
2236            if sol.status != 0:
2237                print(f"⚠️ Integration warning: {sol.message}")
2238            
2239            n_points = sol.y.shape[1]
2240            if n_points < n_frames:
2241                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2242                n_frames = n_points
2243
2244            x_vals, xi_vals = sol.y
2245    
2246            if projection is None:
2247                projection = 'phase'
2248    
2249            fig, ax = plt.subplots()
2250            point, = ax.plot([], [], 'ro')
2251            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2252    
2253            if projection == 'phase':
2254                ax.set_xlabel('x')
2255                ax.set_ylabel(r'$\xi$')
2256                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2257                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2258    
2259                def update(i):
2260                    point.set_data([x_vals[i]], [xi_vals[i]])
2261                    traj.set_data(x_vals[:i+1], xi_vals[:i+1])
2262                    return point, traj
2263    
2264            elif projection == 'position':
2265                ax.set_xlabel('x')
2266                ax.set_ylabel('x')
2267                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2268                ax.set_ylim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2269    
2270                def update(i):
2271                    point.set_data([x_vals[i]], [x_vals[i]])
2272                    traj.set_data(x_vals[:i+1], x_vals[:i+1])
2273                    return point, traj
2274    
2275            elif projection == 'frequency':
2276                ax.set_xlabel(r'$\xi$')
2277                ax.set_ylabel(r'$\xi$')
2278                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2279                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2280    
2281                def update(i):
2282                    point.set_data([xi_vals[i]], [xi_vals[i]])
2283                    traj.set_data(xi_vals[:i+1], xi_vals[:i+1])
2284                    return point, traj
2285    
2286            else:
2287                raise ValueError("Invalid projection mode")
2288    
2289            ax.set_title(f"1D Singularity Flow ({projection})")
2290            ax.grid(True)
2291            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2292            plt.close(fig)
2293            return ani
2294    
2295        elif self.dim == 2:
2296            x, y = self.vars_x
2297            xi, eta = symbols('xi eta', real=True)
2298    
2299            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
2300            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
2301            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
2302            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
2303    
2304            def hamilton(t, Y):
2305                x, y, xi, eta = Y
2306                return [
2307                    dxdt(x, y, xi, eta),
2308                    dydt(x, y, xi, eta),
2309                    dxidt(x, y, xi, eta),
2310                    detadt(x, y, xi, eta)
2311                ]
2312    
2313            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0],
2314                            t_eval=np.linspace(0, tmax, n_frames))
2315
2316            if sol.status != 0:
2317                print(f"⚠️ Integration warning: {sol.message}")
2318            
2319            n_points = sol.y.shape[1]
2320            if n_points < n_frames:
2321                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2322                n_frames = n_points
2323                
2324            x_vals, y_vals, xi_vals, eta_vals = sol.y
2325    
2326            if projection is None:
2327                projection = 'position'
2328    
2329            fig, ax = plt.subplots()
2330            point, = ax.plot([], [], 'ro')
2331            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2332    
2333            if projection == 'position':
2334                ax.set_xlabel('x')
2335                ax.set_ylabel('y')
2336                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2337                ax.set_ylim(np.min(y_vals) - 1, np.max(y_vals) + 1)
2338    
2339                def update(i):
2340                    point.set_data([x_vals[i]], [y_vals[i]])
2341                    traj.set_data(x_vals[:i+1], y_vals[:i+1])
2342                    return point, traj
2343    
2344            elif projection == 'frequency':
2345                ax.set_xlabel(r'$\xi$')
2346                ax.set_ylabel(r'$\eta$')
2347                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2348                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2349    
2350                def update(i):
2351                    point.set_data([xi_vals[i]], [eta_vals[i]])
2352                    traj.set_data(xi_vals[:i+1], eta_vals[:i+1])
2353                    return point, traj
2354    
2355            elif projection == 'phase':
2356                ax.set_xlabel('x')
2357                ax.set_ylabel(r'$\eta$')
2358                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2359                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2360    
2361                def update(i):
2362                    point.set_data([x_vals[i]], [eta_vals[i]])
2363                    traj.set_data(x_vals[:i+1], eta_vals[:i+1])
2364                    return point, traj
2365    
2366            else:
2367                raise ValueError("Invalid projection mode")
2368    
2369            ax.set_title(f"2D Singularity Flow ({projection})")
2370            ax.grid(True)
2371            ax.axis('equal')
2372            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2373            plt.close(fig)
2374            return ani
2375
2376    def interactive_symbol_analysis(pseudo_op,
2377                                    xlim=(-2, 2), ylim=(-2, 2),
2378                                    xi_range=(0.1, 5), eta_range=(-5, 5),
2379                                    density=100):
2380        """
2381        Launch an interactive dashboard for symbol exploration using ipywidgets.
2382    
2383        This function provides a user-friendly interface to visualize various aspects of the pseudo-differential operator's symbol.
2384        It supports multiple visualization modes in both 1D and 2D, including group velocity fields, micro-support estimates,
2385        symplectic vector fields, symbol amplitude/phase, cotangent fiber structure, characteristic sets and Hamiltonian flows.
2386    
2387        Parameters
2388        ----------
2389        pseudo_op : PseudoDifferentialOperator
2390            The pseudo-differential operator whose symbol is to be analyzed interactively.
2391        xlim, ylim : tuple of float
2392            Spatial domain limits along x and y axes respectively.
2393        xi_range, eta_range : tuple
2394            Frequency domain limits along ξ and η axes respectively.
2395        density : int
2396            Number of points per axis used to construct the evaluation grid. Controls resolution.
2397    
2398        Notes
2399        -----
2400        - In 1D mode, sliders control the fixed frequency (ξ₀) and spatial position (x₀).
2401        - In 2D mode, additional sliders control the second frequency component (η₀) and second spatial coordinate (y₀).
2402        - Visualization updates dynamically as parameters are adjusted via sliders or dropdown menus.
2403        - Supported visualization modes:
2404            'Symbol Amplitude'           : |p(x,ξ)| or |p(x,y,ξ,η)|
2405            'Symbol Phase'               : arg(p(x,ξ)) or similar in 2D
2406            'Micro-Support (1/|p|)'      : Reciprocal of symbol magnitude
2407            'Cotangent Fiber'            : Structure of symbol over frequency space at fixed x
2408            'Characteristic Set'         : Zero set approximation {p ≈ 0}
2409            'Characteristic Gradient'    : |∇p(x, ξ)| or |∇p(x₀, y₀, ξ, η)|
2410            'Group Velocity Field'       : ∇_ξ p(x,ξ) or ∇_{ξ,η} p(x,y,ξ,η)
2411            'Symplectic Vector Field'    : (∇_ξ p, -∇_x p) or similar in 2D
2412            'Hamiltonian Flow'           : Trajectories generated by the Hamiltonian vector field
2413    
2414        Raises
2415        ------
2416        NotImplementedError
2417            If the spatial dimension is not 1D or 2D.
2418    
2419        Prints
2420        ------
2421        Interactive matplotlib figures with dynamic updates based on widget inputs.
2422        """
2423        dim = pseudo_op.dim
2424        expr = pseudo_op.expr
2425        vars_x = pseudo_op.vars_x
2426    
2427        mode_selector_1D = Dropdown(
2428            options=[
2429                'Symbol Amplitude',
2430                'Symbol Phase',
2431                'Micro-Support (1/|p|)',
2432                'Cotangent Fiber',
2433                'Characteristic Set',
2434                'Characteristic Gradient',
2435                'Group Velocity Field',
2436                'Symplectic Vector Field',
2437                'Hamiltonian Flow',
2438            ],
2439            value='Symbol Amplitude',
2440            description='Mode:'
2441        )
2442
2443        mode_selector_2D = Dropdown(
2444            options=[
2445                'Symbol Amplitude',
2446                'Symbol Phase',
2447                'Micro-Support (1/|p|)',
2448                'Cotangent Fiber',
2449                'Characteristic Set',
2450                'Characteristic Gradient',
2451                'Symplectic Vector Field',
2452                'Hamiltonian Flow',
2453            ],
2454            value='Symbol Amplitude',
2455            description='Mode:'
2456        )
2457    
2458        x_vals = np.linspace(*xlim, density)
2459        if dim == 2:
2460            y_vals = np.linspace(*ylim, density)
2461    
2462        if dim == 1:
2463            x, = vars_x
2464            xi = symbols('xi', real=True)
2465            grad_func = lambdify((x, xi), diff(expr, xi), 'numpy')
2466            symplectic_func = lambdify((x, xi), [diff(expr, xi), -diff(expr, x)], 'numpy')
2467            symbol_func = lambdify((x, xi), expr, 'numpy')
2468
2469            xi_slider = FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2470            x_slider = FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2471    
2472            def plot_1d(mode, xi0, x0):
2473                X = x_vals[:, None]
2474    
2475                if mode == 'Group Velocity Field':
2476                    V = grad_func(X, xi0)
2477                    plt.quiver(X, V, np.ones_like(V), V, scale=10, width=0.004)
2478                    plt.xlabel('x')
2479                    plt.title(f'Group Velocity Field at ξ={xi0:.2f}')
2480    
2481                elif mode == 'Micro-Support (1/|p|)':
2482                    Z = 1 / (np.abs(symbol_func(X, xi0)) + 1e-10)
2483                    plt.plot(x_vals, Z)
2484                    plt.xlabel('x')
2485                    plt.title(f'Micro-Support (1/|p|) at ξ={xi0:.2f}')
2486    
2487                elif mode == 'Symplectic Vector Field':
2488                    U, V = symplectic_func(X, xi0)
2489                    plt.quiver(X, V, U, V, scale=10, width=0.004)
2490                    plt.xlabel('x')
2491                    plt.title(f'Symplectic Field at ξ={xi0:.2f}')
2492    
2493                elif mode == 'Symbol Amplitude':
2494                    Z = np.abs(symbol_func(X, xi0))
2495                    plt.plot(x_vals, Z)
2496                    plt.xlabel('x')
2497                    plt.title(f'Symbol Amplitude |p(x,ξ)| at ξ={xi0:.2f}')
2498    
2499                elif mode == 'Symbol Phase':
2500                    Z = np.angle(symbol_func(X, xi0))
2501                    plt.plot(x_vals, Z)
2502                    plt.xlabel('x')
2503                    plt.title(f'Symbol Phase arg(p(x,ξ)) at ξ={xi0:.2f}')
2504    
2505                elif mode == 'Cotangent Fiber':
2506                    pseudo_op.visualize_fiber(x_vals, np.linspace(*xi_range, density), x0=x0)
2507    
2508                elif mode == 'Characteristic Set':
2509                    pseudo_op.visualize_characteristic_set(x_vals, np.linspace(*xi_range, density), x0=x0)
2510    
2511                elif mode == 'Characteristic Gradient':
2512                    pseudo_op.visualize_characteristic_gradient(x_vals, np.linspace(*xi_range, density), x0=x0)
2513    
2514                elif mode == 'Hamiltonian Flow':
2515                    pseudo_op.plot_hamiltonian_flow(x0=x0, xi0=xi0)
2516    
2517            # --- Dynamic container for sliders ---
2518            controls_box = VBox([mode_selector_1D, xi_slider, x_slider])
2519            # --- Function to adjust visible sliders based on mode ---
2520            def update_controls(change):
2521                mode = change['new']
2522                # modes that depend only on xi and eta
2523                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)',
2524                            'Group Velocity Field', 'Symplectic Vector Field']:
2525                    controls_box.children = [mode_selector_1D, xi_slider]
2526                # modes that require xi and x
2527                elif mode in ['Hamiltonian Flow']:
2528                    controls_box.children = [mode_selector_1D, xi_slider, x_slider]
2529                # modes that require nothing
2530                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2531                    controls_box.children = [mode_selector_1D]
2532            mode_selector_1D.observe(update_controls, names='value')
2533            update_controls({'new': mode_selector_1D.value}) 
2534            # --- Interactive binding ---
2535            out = interactive_output(plot_1d, {'mode': mode_selector_1D, 'xi0': xi_slider, 'x0': x_slider})
2536            display(VBox([controls_box, out]))
2537
2538        elif dim == 2:
2539            x, y = vars_x
2540            xi, eta = symbols('xi eta', real=True)
2541            symplectic_func = lambdify((x, y, xi, eta), [diff(expr, xi), diff(expr, eta)], 'numpy')
2542            symbol_func = lambdify((x, y, xi, eta), expr, 'numpy')
2543
2544            xi_slider=FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2545            eta_slider=FloatSlider(min=eta_range[0], max=eta_range[1], step=0.1, value=1.0, description='η₀')
2546            x_slider=FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2547            y_slider=FloatSlider(min=ylim[0], max=ylim[1], step=0.1, value=0.0, description='y₀')
2548    
2549            def plot_2d(mode, xi0, eta0, x0, y0):
2550                X, Y = np.meshgrid(x_vals, y_vals, indexing='ij')
2551    
2552                if mode == 'Micro-Support (1/|p|)':
2553                    Z = 1 / (np.abs(symbol_func(X, Y, xi0, eta0)) + 1e-10)
2554                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='inferno')
2555                    plt.colorbar(label='1/|p|')
2556                    plt.xlabel('x')
2557                    plt.ylabel('y')
2558                    plt.title(f'Micro-Support at ξ={xi0:.2f}, η={eta0:.2f}')
2559    
2560                elif mode == 'Symplectic Vector Field':
2561                    U, V = symplectic_func(X, Y, xi0, eta0)
2562                    plt.quiver(X, Y, U, V, scale=10, width=0.004)
2563                    plt.xlabel('x')
2564                    plt.ylabel('y')
2565                    plt.title(f'Symplectic Field at ξ={xi0:.2f}, η={eta0:.2f}')
2566    
2567                elif mode == 'Symbol Amplitude':
2568                    Z = np.abs(symbol_func(X, Y, xi0, eta0))
2569                    plt.pcolormesh(X, Y, Z, shading='auto')
2570                    plt.colorbar(label='|p(x,y,ξ,η)|')
2571                    plt.xlabel('x')
2572                    plt.ylabel('y')
2573                    plt.title(f'Symbol Amplitude at ξ={xi0:.2f}, η={eta0:.2f}')
2574    
2575                elif mode == 'Symbol Phase':
2576                    Z = np.angle(symbol_func(X, Y, xi0, eta0))
2577                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='twilight')
2578                    plt.colorbar(label='arg(p)')
2579                    plt.xlabel('x')
2580                    plt.ylabel('y')
2581                    plt.title(f'Symbol Phase at ξ={xi0:.2f}, η={eta0:.2f}')
2582    
2583                elif mode == 'Cotangent Fiber':
2584                    pseudo_op.visualize_fiber(np.linspace(*xi_range, density), np.linspace(*eta_range, density),
2585                                              x0=x0, y0=y0)
2586    
2587                elif mode == 'Characteristic Set':
2588                    pseudo_op.visualize_characteristic_set(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2589                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2590    
2591                elif mode == 'Characteristic Gradient':
2592                    pseudo_op.visualize_characteristic_gradient(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2593                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2594    
2595                elif mode == 'Hamiltonian Flow':
2596                    pseudo_op.plot_hamiltonian_flow(x0=x0, y0=y0, xi0=xi0, eta0=eta0)
2597                    
2598            # --- Dynamic container for sliders ---
2599            controls_box = VBox([mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider])
2600            # --- Function to adjust visible sliders based on mode ---
2601            def update_controls(change):
2602                mode = change['new']
2603                # modes that depend only on xi
2604                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)', 'Symplectic Vector Field']:
2605                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider]
2606                # modes that require xi, eta, x and y
2607                elif mode in ['Hamiltonian Flow']:
2608                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider]
2609                # modes that require x and y
2610                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2611                    controls_box.children = [mode_selector_2D, x_slider, y_slider]
2612            mode_selector_2D.observe(update_controls, names='value')
2613            update_controls({'new': mode_selector_2D.value}) 
2614            # --- Interactive binding ---
2615            out = interactive_output(plot_2d, {'mode': mode_selector_2D, 'xi0': xi_slider, 'eta0': eta_slider, 'x0': x_slider, 'y0': y_slider})
2616            display(VBox([controls_box, out]))